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ABSTRACT 


The problem of optimal load flow in power systems is 
considered. Based on Carpentier's formulation, the problem is re- 
formulated to achieve a reduction in the number of system variables, 
by treating all generator nodes as swing busses. Reduced system 
variables are also classified into control and dependent variables. 
The set of equality constraints is associated with the objective 
function using Lagrangian multipliers, and functional inequality 


constraints are included as penalties. 


The set of control variables is iteratively adjusted, 
using Newton's method, to minimize the objective function, while 
the set of dependent variables is evaluated by solving the set of 


equality constraints after each such adjustment. 


The reduction in the number of system variables and their 
classification as such provide a great saving in computer storage 
requirement as compared to other established methods. Furthermore, 


the use of Newton's method provides an excellent convergence behaviour. 


The developed minimization algorithm is applied to the 
well knovn minimum generation cost and minimum system loss problem. 
In addition, the minimum fuel consumption, and combined fuei-cost 


minimization problems are defined and solved. 


iv 


Asai voces 2} eoidatrey fongnos ¥6 $98 ont | 
AttosAd axtntntn of ybodtem 7 netwelt gated 


‘ighin andar danes-oaee mt wotdouber SAT 


sesose: vadugnos at’ paige Seung s Shivong Wout es rotYeotthazela 


An optimal ordering scheme of system nodes, for use with 
large systems, is also developed and compared to two other effective 
schemes. The new scheme proved to be generally comparable to both 


schemes for the cases studied. 


The concept of fixed penalty factor, developed in assocation 
with the minimization algorithm is investigated and compared to 
the usual concept of monotonically increasing sequence of penalty 
factors. The former proved to be superior from a practical point of 


view. 
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CHAPTER I 


INTRODUCTION 
With the rapid growth of demand of electrical energy, 
power systems have had to expand. Associated with such expansion 
are problems which have not been encountered before, or encountered 
only to a minor extent. These problems are related to the question 
of how to operate or expand a power system reliably and economically. 
Today, with the vast sizes of power systems that exist, many, if not 


all, of these problems still have to be faced and dealt with. 


While some of these problems, such as stability and associated 
problems, are related to the transient behaviour of the power system, 
some others are related to the steady state operation of the system. 
These latter problems include, among others, the problem of optimal 
load flow, i.e. the problem of obtaining a load-flow solution which 


optimizes a certain operating criterion. 


The importance of the optimal load flow problem is evident 
if one considers the very simple case of two generators supplying a 
common load, and recognizes that cost savings can be effected if 
generation is shifted from the less efficient generator to the more 
efficient one. However, an immediate question arises: how much generation 
is to be shifted to accomplish the most savings? This is the heart of 


the optimal load flow problem. 
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1.1 Development of the Optimal Load-Flow Problem 


The optimal load flow problem was first recognized in a 
restricted sense in the form of the economic real power dispatch problem. 
As mentioned before, the aim was to obtain a generation schedule that 
will meet the demand most economically. Generator voltages were kept 


fixed, thus no control on reactive generation could be effected. 


First attempts to solve this problem neglected the system's 
transmission losses. The result was a generation schedule that requires 


all generators to operate at the same incremental égxbl, 


The advent of integrated systems and the interconnection 
between power companies for the purpose of economy interchange, rendered 
this approach unsatisfactory. Transmission losses became significant, 
and their effect had to be taken into account. This led to the 
modification of the incremental cost method using penalty factors and 


loss Pea tia debt aGeual’ 


Although exact coordination of incremental loss and incremental 
cost was AER IeVEdE bee the method suffers from the fact that the methods 
used to determine these loss coefficients (B - constants) are approximate 


in nature and depend on particular load levels and generator voltages. 


Relatively recently an exact transmission loss formula was 
developed using the bus impedance matrix of the system, and the resuits 
of a load flow ee titionens-s Transmission loss was expressed in terms 


of real and reactive powers at all nodes. Coefficients of the formula 
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are not constants but change as node voltages and angles change. The 
formula was used in conjunction with the exact co-ordination equations 


to produce the exact solution of the economic real power dispatch problem. 


In the same reference, an extension to the method to 
economically allocate reactive generation is presented. By moving along 
the negative of the gradient of transmission loss with respect to 
reactive generations, these losses can be reduced. By alternating real 
and reactive power dispatch processes, it is claimed that the more 
general problem of exact economic dispatch can be solved. However, it 
has been later shown that although,this extension, will allow the 
control of reactive generation, it will fail to locate the optimum 


[13], 


solution 


1.2 The General Optimal Load Flow Problem: Carpentier Formulation 


In 1962 a breakthrough in the problem of optimal load flow 
was obtained by parpeniicree 2, who developed an exact and general 


formulation of the economic dispatch problem. 


An N bus power system is considered. At each node there may 
be generation Py and Qg> and consumption C and D, of active and reactive 
power, respectively. Each node is also characterized by its voltage V, 
and phase angle 6. These variables are interrelated, in a nonlinear 
fashion, by the network relations which govern the flow of power in 


every and all parts of the system. 


Of course, total production ) P 


and ) Q. , must equal the 
i yl a 
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total consumption ) C. and } D., plus the respective transmission losses. 
j i 
This can be guaranteed by the satisfaction of the above mentioned network 


relations. 


Further, due to engineering and physical considerations, 
some variables can only change within a specified range. For example, 
a generator can not produce any power beyond a certain value determined 


by the capacity of the boiler and turbine driving it. 


Finally, the cost of operating the power system is a function 
of real power generations Py » but not of reactive generations Q : 
i j 
which are cost free once the equipment required for their production has 


been installed. 


Thus the problem consists of minimizing the operating cost 
»--e5P ) Subject to the following: 
9 aN 
1) Satisfaction of equality constraints formed by the power 
flow equations (network relations). 


P. =C. = , iF Mo Vito cos(s.-6 -8,) (eta ) 


58; 
where Y.e ie is the term in the node admittance matrix 


corresponding to nodes i and a. 
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2) Satisfaction of inequality constraints imposed by the 


Operating limits of the various variables: 


Sh Ds (Micador (13) 
Thin 95 Tmax 
Es iid ts pees 101 (1.4) 
Thin g; thax 
V. Es ee Ges) 
‘nin ! Tmax 


Other inequality constraints can also be imposed. 


Although Carpentier presented the problem in the form of 
economic dispatch, changing the function f to any other power system 
function (e.g. total real generation), and using the appropriate 
constraints, will result in a different power system optimization 


problem of the same general form. 
1.3 Optimal Load-Flow Solution Methods 


Apart from the methods that solve only special cases of the 
optimal load flow problem, some of which have been described in Section 


1.1, several methods have been developed to tackle the problem in its 


most general ‘anglaise A comparison between some of these methods 


can be found in reference 13. 


In this section two of these methods, due to Dommel and 


Tinney eee and Sasson, Viloria and Aboytest 221 , will be described in 


some detail because of their relevance to the work presented in this 
thesis. However, a brief outline of the other methods will be given first. 


ede, Alt6 517.) and Ramamoorty and Raol 18] employed nonlinear 


programming methods to solve the problem. Several techniques have been 


described for incorporating the problem constraints as penalties on 
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the cost function, thus transforming the problem of constrained 


optimal load flow to an unconstrained problem. 


While in reference 18 a first order gradient technique was 
used for the unconstrained minimization process, references 16 and 17 
used the Fletcher-Powel ] algorithm 224 which represents cne of the 
most powerful techniques in what is known as the “variable metric 


rathode co 


Though this method worked quite well for small systems, 
it developed convergence and computer storage problems as system size 
increased, and decomposition techniques had to be usedto>d, 


Billinton and Sachdeval !9+201 


used a suboptimal technique 

for solving the problem. With an assumed real power schedule, and 
keeping generator powers and relative phase angles between system nodes 
fixed, voltage magnitudes are optimized using transmission losses, a 
function of the voltages only, as an objective function, reactive power 
equations as equality constraints, ard the limits on reactive generation 
and voltage magnitudes as inequality constraints. With the resulting | 
voltage magnitudes fixed, a real power dispatch is obtained optimizing 


the cost of generation under the usual equality and inequality constraints. 


The process is then repeated until no further improvements can be obtained. 


Although the method is similar to that of reference 12, the 
way the problem is formulated guarantees an optimal solution. However, 


[20] 


this optimal solution will depend on the choice of the swing bus 


Reid and Hasdortfee. applied quadratic programming to solve 


stata okey" iain debian iutrewog see ) 

enaiaye Flame vol (Tew atop: hettnow Waryon aFAr duet OS) ypodiom : 

site made as aan suceie Sato boo vonagravnes beqataveb sh 

“od bot. auphmtom satst eoqiuosd on sbeanstaa 

 sayphtsna ther scala baa (080F I, atte bins andthe va a 

bis safubrloe veWoq (Phen omyiae ut; NM Lonidove 9 gntv6a ON 

debon meteye nsowibd calenaieandn oytzetin bent zron0g woheraday parietal - 

8 upeatal abl eatnenets ‘prion patna: ss abut ng aossioy bax 

‘aWoq, vinden «oF tou avbaowtido ah Zt Xin expat ivy ont to noi sant. 

edt sevens 2-500" SHIM 942 ow <2dubgea eno pre 

won saa 9 ail AFM, «antheWeives (attained 2h exbud 

er biiniana! Wont eo 2¥ dataaethy show Toon 6° ith eebutt gem 

Leanisrdenno Ari acpant bie rie (ouav ait “ati nor aorenae to-2209 8 

maaescandbaler sia hc pidial in ee ed 2 % 
a Lane 7 “is ma / ; 

“ fare comers a. — a hike ‘ 


Cue 
7 - 


, 


7 ae i Bays 


the problem. Since the method requires a quadratic cost function and 
linear equality constraints, new variables had to be introduced to 
transform system operating cost into a quadratic function of system 
variables including those newly defined. These new variables also 
includes the variables introduced to transform the inequality 


constraints into equality constraints. 


The method worked very nicely for systems of up to 118 
busses in size. No mention of storage requirements was given, however. 
But since the number of variables involved far exceeds those of any 
other method, storage requirements can jeoparadize the method's success 
for large realistic systems. 


gules centres around ordinary load 


Dommel and Tinney's metho 
flow solution by Newton's mephon ooo Lagrangian multipliers are used 
to associate the equality constraints with the objective function. 


Inequality constraints are included as penalties. 


The equality constraints consist of all equations forming 
the ordinary load-flow problem, i.e. one real power equation for each 
generator node, and two real and reactive power equations for each 


load node. Slack node equations are not included. 


Apart from fixed parameters and those which can be readily 
obtained from given equations, e.g. reactive generations, system 
variables are real power generations, with that of the slack node 
expressed in terms of the system voltages and angles which form the rest 
of the variables. The variables are classified as control or 
independent variables "u" consisting of generator voltages and real 


powers, and dependent variables "x" formed by voltages and angles of load 
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nodes as well as angles of generator nodes. 
The Lagrangian function is thus formed as: 
L(usx) = f(usx) + a! g(u,x) (1.6) 


where f(u,x), g(u,x) and A are the objective function, equality 
constraints and the Lagrangian multipliers, respectively. Conditions 


for the minimum of the objective function are: 


oL 

x16 g(u,x) = 0 (iar) 
ag 

aL at a belie 

ax ax * [ sx 1 RS= 0: (1.8) 
og 

ab _ af agi 

le ee eo a (1.9) 


Given an estimate for the control variables "u", equations 
(1.7) which are the load-flow equations, are first solved for the 
dependent variables "x". The Jacobian is then used in (1.8) to solve 
for "A". Equations (1.9) will then give the gradient VE of the 
Lagrangian function with respect to the control variables u. A 
correction is then applied to these control variables by a move along 
the negative direction of VE which is the direction of maximum decrease 
of Lat that point. This correction is given by: 


ni Se oon (1210) 


where c is an acceleration factor determining the size of the move. 


The process is repeated until no further improvements can be 
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The method, although theoretically sound, suffers from the 
poor convergence characteristics of the steepest descent method and 
its sensitivity to accelerationc. Storage requirements could also 
pose a problem for large systems, due to the large number of equality 
constraints, even if only non-zero elements of the Jacobian matrix are 
stored. Another problem is that the solutions obtained for some 
problems, e.g. the minimum loss problem (f = P (V,6)), are dependent 


slack 
on the location of the slack node. 


The major drawback of the method is, however, its inability 


[27] 


to handle complex objective functions » possibiy because of the 
inadequacy of the steepest descent direction. Complexity of the cost 
function can result from penalty terms introduced by violated nonlinear 
inequality constraints which badly deform the hypercontours in the 


state space. This is possibly why no results involving such constraints 


have been published using this method. 


sassont““daid not distinguish between equality and inequality 
constraints, and penalized the objective function for both types. 
Generator equations were first removed fron the equality constraints and 
substituted into the objective function and inequality constraints to 
eliminate generator real and reactive powers as variables. To minimize 
the cost function, corrections to the voltages and angles of all nodes 
are applied at the same time using the Hessian matrix of second order 
partial derivatives of the panalized objective function. 


au=-H'yv f (1.11) 
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where u, H and an are the vector of voltages and angles, the Hessian 


matrix and the gradient of the penalized objective function with respect 


to the vector u, respectively. 


The method requires even more storage than Dommel and Tinney's 
method, thus its success could be limited for large systems. Since 
equality constraints are not satisfied until the last iteration, they 
are always present as penalties in the objective function with the 
obvious adverse effects on the method's convergence. This can be seen 
from the jess than adequate satisfaction of these equality constraints 


after a solution is claimed to be obtainedeoc4. 
1.4 Research Objective 


As shown in Sections 1.1 and 1.3 every method suffers from 
one drawback or another. This limits their use in the power industry 
either due to their storage requirements or poor convergence. Thus it 
was considered worthwhile to investigate the possibility of developing 
another method which possesses the advantages of the previous methods, 


and at the same time does not have, at least, their major shortcommings. 


Also investigated, is a new approach to optimal ordering of 
system nodes to minimize the number of new off-diagonal elements 
introduced during the elimination process used in Newton's method. 
The idea is to locate as much off-diagonal element of the original 


matrix as possible in the lower right hand corner of such matrix. 


The goal of this work is, thus, to develop a method and its 
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Supporting mechanisms that can be readily implemented by the power 


industry. 
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CHAPTER II 


PROBLEM REFORMULATION AND ITS SOLUTION 

As mentioned in Section 1.4, the objective of the solution 
method is two fold; first, to reduce computer eres requirement, 
and second, to obtain better convergence behaviour. The first objective 
can be achieved by reformulating the problem to eliminate some of the 
variables, and using a solution method that does not handle all the 
remaining variables simultaniously. The second objective can be 
fulfilled by choosing a minimization criterion which provides a better 


minimization direction. 
2.1 Problem Reformation 


The basic Carpentier formulation was stated in Section 1.2. 
It was also mentioned that although Carpentier presented his formulation 
in the form of the economic dispatch problem, changing the objective 
function will lead to a different power system optimization problem. 
This means that the objective function f is not bound to be a function of 
real power generations only. Therefore, it would be preferrable to 
restate the problem in the following form, which differs from that of 
Section 1.2.in only the definition of the objective function: 

Minimize the scalar function f of system variables subject 


to the following constraints: 
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Equality constraints: 


P. (V,6) 3 Veusigy aa’ (2.1) 
aie” Naseem pall 
Q. (V6) ii (0,93) . 4s (2pi2,) 
Inequality constraints: 
P Sols egy (23) 
“hin 9; max 
i for generator 
busses, j for all 
Q. 20s 0. busses (224) 
min 9; max 
Ve Eon corcan | (255) 
Jinin J Jmax ) 


where, at a bus i, P. (V8) and Q. (V6) are reali and reactive power 
injections given by the right hand sides of equations (1.1) and (1.2) 


respectively. 


The equality constraints (load flow equations) consist of 
2N equations. This number can be reduced if those equations corresponding 
- to generator nodes are removed. This can be accomplished by substituting 
for ee and aoe of generator nodes in the objective function and/or 
the appropriate inequality constraints. This leaves two equations per 
each load node to form the equality constraints. This is gh le 
equations less than an ordinary load flow equations set, where Ng is 
the number of generator busses. Note that in an ordinary load flow 


problem, reactive power equations of generator busses are always excluded 


due to fixed generator voltages. Also excluded is the real power 
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equation of the slack node due to the fact that its angle is fixed as 


reference, and its power is left floating to absorb system losses. 


Elimination of PG and Qg in such a way, and the fact that 
j 7 
almost all system variables can be expressed in terms of system voltages 


and angles, transform the problem into the following form. 


Minimize the scalar function f(V,6) subject to the following 
constraints: 
Equality constraints at a load bus 7: 


Ejeet Ge = 0 (2,6)% 
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Inequality constraint at any bus k: 


VSM SM (2.10) 


Kin max 


* Note that Ea and Qa are zero at a load bus. 
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The prime on the upper and lower limits of generator powers 
indicates that these limits are modified to account for local loads. 
It will be omitted from now on,but the foregoing should be emphasized. 
Note also that the main variables of the system have been reduced to 


include only the voltage magnitudes and angles. 


Elimination of 06 and ta also means that all generator 
busses are treated as swing busses in the sense that their powers are 
determined from the voltage distribution in the system. This will 
eliminate the dependence of the solution of some problems, e.g. the 
minimum loss problem, on the choice of a particular swing bus. Any 
One of these busses can be chosen as a reference bus and its angle set 
to zero. Note here tnat the voltage magnitude at any of these busses 
is free to change (within prescribed limits) and the voltage level in 


the system is no longer determined by a fixed voltage at a swing bus. 


Load flow equations (2.6) and (2.7) are, in vector form: 

g(V,6) = 0 (2a) 
Thevorder ony themvectortg 19"2 Nis where NPEN-N, is the number of load 
busses. Solution of (2.11) will provide the values of en, unknowns out 
of 2N-1 variables, thus the remaining 2N I variables should be assumed. 
This provides the basis of variable classification into: 

a) Control or specified variables; 2Ngo in number, 

b) Dependent or unknown variables; en, in number. 

The most logical and natural choice is that the first set, denoted by the 


vector u, should include the voltage magnitude and angles at generator 
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busses. They are exactly aNaol in number (note that the angle of the 
reference bus is fixed at zero). The second set, denoted by the vector 
x, will consist of the voltage magnitudes and angles at load busses. 
The number of these quantities is exactly eN 

Thus, load flow equations (2.11) become 

g(u,x) = 0 aes Wa 
and the statement of the problem takes the form: 

Minimize the scaiar function f(u,x) subject to: 

Equality constraints: 

g (u,x) = 0 (2120 


Inequality constraints; 


ce ek rn (2513) 
min max 

ree Dea Us X Jee (2.104) 
DIC etiam ed ) | 

V.. 2a Ee arm 215 
Jmin~ 37 Jax 


where i iS a generator bus, and j is any bus. 


2.2 Solution of the Optimal Load Flow Problem 


Any constrained nonlinear programming problem, like the one 
described in the previous section, can be solved by converting it into an 
unconstrained problem, and then applying one of the numerous methods 
developed for unconstrained minimization. The conversion is achieved 
by defining an appropriate auxiliary function, in terms of the original 


problem functions, and using it as a new unconstrained objective function. 


There are two ways to incorporate the problem constraints 
into the objective function to form this new auxiliary function. The 


first is to use the Lagrangian multiplier theorem for the equality 


at ‘te alana git Jand ode ni hy 3s vite, oH eo a a 


— wetkaw ott vd by tort .tae bidase itt pak dn box ot tab sone 
ad \ 

toeeud bwol 36 nares hie gobuytinpan spetfov ary vo satanen, hi 
“gis yhyases 21 earatidaup oad . ai mun C 


essed (TPS) 2enolhtsups «af? ‘baal euetT : 


( = “ ) ; as eal 


“ , : - 

‘wat Sat eon asfdoug any Yo Sen siod2 only 
7 
j ! 


(apn)? poryanu? Wisse wm auto aii 


-stiatewewo yotioupa 


($t.9) ; o' xo) 
-ejoteiiens yattewpeal . 
AB) pe wt ae 
ory 
red was ey D m4 
(ar) 4M » > :¥ 


|| ae a vein 
,eud vas ef & be , eee Oso PNOD B et te 


mabdir’ volt bec. Tembban og te sila . 


uno aft sdtl ,de lider onto Year laon boatevtenes ww 8. 

ne obnt 4% pobtedvann wi bavloe ad ama. .aoriose 2yolyvearg any at bed! tpaes 
4 

sbonsan eporsmyn att To snd gad age melt bus nat oote cantantt 


bavetion 2! noterayaos oT .nolsastafarn bentaWenognn der 


2 


_ Teataine sid to eirist al .notioa we ierk raves pean 
uesineeianieniciete worn 26 SF hi ee 200 sano ms 
. 


‘ ene 4 a open Lone ions 
bm eh 0 148 ot 


a ? 


say * OT its ogy vata netg yerp on in 


r 


aie 


17 


constraints, and the Kuhn - Tucker syenen oe for inequality 
constraints. This requires that the equality constraints are explicitly 
Satisfied. Also required is the satisfaction of what is known as the 
"exclusion equations" of the Kuhn-Tucker theorem- '31) The minimization 
process thus moves from one feasible point to another until the 


minimum of the objective function is located. 


The second way is to penalize the objective function for both 
types of constraints, as Sasson et. al. didt22], Many forms for the 


[24,28,29] In this 


penalty function are available in the literature 
case, one moves in the whole space rather than the feasible region until 


the optimum is obtained. 


To choose between these two approaches, one should recognize 
that unless the equality constraints are explicitly satisfied, they will 
always be present as penalties with the obvious adverse effects on 
convergence as mentioned before. Thus the idea of penalizing the cost 
function for this type of constraint is ruled out. Furthermore, very 
few inequality constraints are violated simultaneously, so their 
inclusion as penalties is not as bad as it would be in the case of 
equality constraints. Moreover, the evaluation of Kuhn-Tucker dual 
variables and the associated change of node type, which would increase 
the number of equality constraints, and the order of the load-flow 


portion, is avoided. 
For the problem at hand, the Lagrangian function is defined 
as: 


L(u,x) = f(u,x) + A g(u,x) (2.16) 
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where, f(u,x) is the objective function which incorporates the penalty 
terms arising from the violated inequality constraints (if any) as will 
be discussed later, and A is the vector of the Lagrangian multipliers 


associated with the equality constraints. 


It is imperative to say that minimization of the Lagrangian 
LZ with the equality constraints satisfied means in fact the minimization 
of f. Thus at the minimum of f(u,x), the following conditions must be 


satisfied giving the optimal solution: 


so = g (usx) = 0 (2.17) 
ag 
ab= ete Taso (2.18) 
ax oS as 
9g 
Sy my Uh eee Se Mea HP | 
dU au [ su J A. 0 (2519) 


Although equations (2.16) - (2.19) are identical to equations 
(1.6) - (1.9) of Dommel and aes one should emphasize the 
difference in the order of matrices and vectors involved, as well as 
the differences in the definition of these various quantities. This is 
what gives this method the advantage as far as computer storage is 


concerned as will be detailed later. 


Equations (2.17) - (2.19) are nonlinear and iterative methods 


are necessary to solve them. The generalized Newton's method has been 


known to be the most powerful of all minimization techniques’?2, resulting 


in a superior convergence behaviour. The only theoretical difficulties 
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associated with this method is the amount of storage required for the 
Jacobian matrix and the expensive efforts to invert it. However, these 
difficulties are drastically reduced in power system studies because 

of the extreme sparsity of the Jacobian so that Newton's method has 


become a standard procedure in the power industry. 


Thus, based on this method, the solution algorithm is as 


follows. A flow chart is also given in Fig. 2.1. 


1) An arbitrary set of values is assumed for system 
voltages and angles. 
2) Load flow equations (2.17) are solved using Newton's 


methodt261 | 


The order of this problem is 2N 
3) Lagrangian multipliers A are obtained from equation 
(2.18) as: 


ag. bas 
Sp | Dead Thoma A (2.20) 


3g 
The equations are linear in d and the Jacobian matrix [~] 
is already available from step 2). If the Jacobian is 


available in factored form (upper and lower triangles)* 


“This is computationally equivalent to the inverse or transposed 


[30] 
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Figure 2.1 Flow-Chart of Solution Algorithm 
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as is usually the case for large systems, this step 
represents one repeat solution of a linear system. 

At this point the values of u, x and » satisfy equations 
(2) )/)eand)(2-18)\ean (hey will) not.wineceneral. 

satisfy equations (2.19) which will, then, give the 
gradient Vf, of the objective function f with respect 
to the eae) variables u when equality constraints are 
satisfiedt!9J, Thus a correction Au in u is possible 
using the relation: 

Au = - nel 


aft (RA) 
where H is the Hessian matrix of the second order partial 
derivatives of the Lagrangian function with respect to 
Usewlts order 1s eNg-1: Formula (2.21) is identical 

to Newton's formula. It gives not only the direction of 


the move, but also its size. 


The new values of the control variables u is then given by 


sep usainie (2.22) 
where the superscript indicates the iteration number. 


If some convergence criterion is satisfied, the solution 


has been found otherwise return to step 2). 


For convergence criterion, one may use one of the following: 


a) 


The change in the value of the objective function is 
less than a prescribed value. The accuracy of this 


criterion is doubtful in cases where the objective function 
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happens to be very flat near its minimum resulting 
in a suboptimal solution. 

b) Each component of the correction vector Au is smaller 
than a preassumed value. This will produce the 
most accurate results because one is checking the 
movement of the operating point itself. 

c) Each component of the gradient vector Mi is less 
than a set tolerance. This is the yee equivalent 
to the theoretical requirement of zero gradient which 
is the condition for an optimal solution. However due 
to the premultiplication by the inverse Hessian and 
the fact that, in practice, the gradient will never be 
identically zero, the correction vector Au may not 


satisfy criterion (b) if it were used instead. Moreover, 


care must be taken because a component Ta of the 
i 


gradient vector Vf will not be zero i Ee RCONLG | 


[15]. 


variable VF 1SHONe ONG) Olga tS aL ts 
2.221 Solution of Equation (2.17) 


Newton's method is well known and its development into the 
most powerful methods of load flow solution by exploiting the sparsity 
of the Jacobian matrix of load flow equations is well cocmmentedsac 
The fact that the Jacobian matrix should be positive definite for assured 


convergence seems to be satisfied for practical power systems as there 


is no evidence to the contrary although very many different systems 
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have been studied. 


In equations (2.17) the vector u is fixed at the initial 
guess, or the corrected value. The vector x is set at an initial 
estimate, or the value from the previous iteration, and the vector 
g(u,x), which is known as the power mismatch, is computed. Generally, 
it will not equal to zero. Corrections in the vector x are then 


carried out using the relations 
Aen SEIU MIG (Oh x!) (2.23) 


SN an rer (2.24) 


th iteration, and 


where a is the Jacobian matrix evaluated at the i 
the superscripts indicate the iteration number. The process is 
repeated until the mismatches g are less than a prescribed value 


(e.g. 107°). 


Equations (2.23) are written in more detail as: 
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where all subscripts belong to the set of load nodes, and P; and q; 
are real and reactive mismatches (forming the vector g) at load node 


i. The elements of the Jacobian matrix are given in Appendix A. 


Even though the Jacobian matrix of load flow equations is 
extremely sparse (3.5% full for the IEEE 118-bus test system). due 
to the fact that each node is connected to only few adjacent nodes, 
the computer time required to invert it in the case of a large power 
system would be prohibitive. Moreover, the inverse proper will be a 


full matrix requiring a large computer storage area. 


These problems have been alleviated by the use of Gaussian 
elimination resulting in an upper triangular matrix, then using 
back substitution to obtain the solution Ax. Forward operations may 
be stored in the lower triangle for the purpose of repeat anne ewe 
It has been shown that for an nxn full matrix the number of operations 
(multiplication-addition) required for triangulization is of the 
order eo compared to n? for proper Tnvered one Back substitution 
would require the same number of operations, — as the premultiplication 


with the inverse. So it is evident that large savings in computer 


time can be effected. 


*=see closure of ref. 2c. 
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25 
This approach combined with an optimal ordering scheme of 
the nodes to preserve the sparsity of the Jacobian during triangulization, 
and the use of compact storage techniques of sparse matrices, in 
which only non-zero elements and their positions are stored, will 
also result in great savings of computer storage which would be 


required otherwise. 


MrosesOVOCiOn Of baudcions (2010) 


After equations (2.17) have been solved, a by product is 
the Jacobian matrix. The vector oF can also be computed. Its 
components are given in Appendix Naas, de cCOSC TUNG LONE WiC mts 
quadratic in real power generations. The availability of these 
two quantities means that, at the present operating point, equations 
(2.18) are linear in the Lagrangian multipliers ». They can be 
solved by virtue of equation (2.20). Whether the Jacobian is available 
in factored form or as explicit inverse Jacobian, the determination 
of A amounts to a repeat solution of the transposed system. In this 
case the vector of replaces the mismatch vector g. The number of 
operations pecan in either case is né if UNee Ma citi xed Sa Ud We coNg 


much lower if the sparsity of the Jacobian is taken advantage of, 


i.e. the factored Jacobian is used instead of the inverse. 


2.2.3 Ihe Hessian Matrix 


In equations (2.21) the Hessian matrix of the second order 
partial derivatives of the objective function with respect to the 
control variables u, is used to compute the corrections Au. The matrix 
is symmetrical, and, thus, only the diagonal and upper triangle 


elements need to be stored. Moreover the matrix is extremely sparse 
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as shown in Figure 2.2 This is due to the fact that the control 
variables consist of the voltages and angles of generator nodes only. 
Unless there is a direct tie between two such nodes, the Lagrangian 
function will be free from any term involving the variables at these 


nodes simultaneously. Hence, no cross term (off diagonal), e.g. 


Z 
oor. will appear in the Hessian except in the case where there 
Lee 
is a tie between generator i and generator j. In practice these ties 


are rare, and the Hessian reduces almost to 2 x 2 submatrices, along 


Na L 3 iB 3 L 
the main diagonal, representing the terms ae Seay 5 See 
5 36. jedi aV 536, 


and ot In Figure 2.2 node j is connected to node k whereas all 
oV 
other nodes are free from any connection with other generator nodes. 


Furthermore, the elements of the Hessian matrix are very 
simple and easy to compute (see Appendix A). This, together with the 
previous analysis of its structure, shows that equations (2.21) are 
easy to handle, and no problem as far as storage and time requirements, 


has to be faced. 


As mentioned earlier, equation (2.21) is identical to 
Newton's relation (2.23). Thus,one expects a convergence behaviour 
as good as Newton's method. The Hessian matrix should be positive 
definite to assure such convergence, a condition which was true in all 
cases studied. The positive definiteness of the Hessian can also be 
assured if the starting point is close to the optimal solution. This 


requirement is always satisfied in a flat starting point*. The fact is 


a 


* A flat start means all voltages are set to its specified value or at 


1.0 p.u. and all angles set at zero. 
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that in any power system the voltage magnitude at any bus will never 
change beyond the range 0.9 - 1.1 p.u., and the phase angles amount 


to a few degrees away from the zero reference. 


2.3 Inequality Constraints 


In the formulation of section 2.1, system variables were 
classified into control and ised variables. As their name 
indicates, the value of any control variable is under control by 
virtue of equations (2.21) and (2.22) or any other valid regulating 
condition. On the other hand, there is no direct control over the 
dependent variables, as they are completely determined by the load 
flow solution. Also, control of the value that any generator power 
can assume is eliminated by treating all generators as swing generators. 
Therefore, inequality constraints, which determine the operating range 
of all these variables, can be classified into parameter (or control) 


and functional constraints. 
2.3.1 Parameter or Control Constraints 


These are the linear inequality constraints applicable to 
each individual control variable. Typical constraints are upper and 


lower limits on voltage magnitudes at generator nodes. 


ihisstype of constraints can be easily satlistiedsby vestrict— 
ing each variable to its operating range, i.e. no control variable 


is allowed to exceed its limit. Thus equations (2.22) are replaced 


by: 
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j j 


ue ite ete AU as Une 
Jnax J J Jmax 
yitl = u ae ul! + Au! U (2527) 
. = a 6 nr kK a ° 
: Jinin J J Snin 
j j ; 
Ula am Nic 
j = otherwise 


The use of this approach not only will keep the successive 
moves within the feasible region as far as control variables are 
concerned, but will also reduce the number of inequality constraints 
to be treated by the penalty function approach, thus reducing the 


chance of introducing penalty terms into the objective function. 
2.3.2 Functional Constraints 


These consist of all problem constraints except parameter 
constraints. In general they are the inequality constraints, linear 
or nonlinear, involving dependent variables and/or two or more of the 
control variables. Typical constraints are upper and lower limits on 
generator powers, and upper and lower limits on voltage magnitude at 


load nodes. 


As mentioned before, the penalty function approacn is used 
to handle this type of constraints,particularly because very few 
such constraints are simultaneously violated. This, in itself is a 
reason to rule out the use of interior point penalty Reon 
which requires that all inequality constraints, whether satisfied or 


not, be monitored and their derivatives calculated. Another reason to 


back this decision is the requirement of a feasible point to start an 
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interior penalty function problem. To locate such a point for a 


high order system is indeed a very tedious task. 


This leaves the exterior point merhodelc 24 as the suitable 
way to handle these functional constraints. For one, they do not 
need a feasible point to start the minimization process. Secondly, 
the derivatives of only the violated constraints need to be computed. 
Finally, in the present practical problem, as may be the case in most 
practical problems, a functional constraint is seldom a rigid limit 
in the strict mathematical sense but is, rather, a soft limit. For 
instance, V<1.0 p.u. on a load bus means that V should not exceed 
Te OsDsuUse bY etoo much. wand, V=I-70] "may still be permissible. Exterior 


point methods do just that (see Figure 2.3) 


Rigid Limits 


Soft Limits 
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Xmin Xmax 


Figure 2,3 Exterior Penalty Function 


To handle inequality constraints (2.13) - (2.15), each 
such relation should be split into two inequality relations, one of 


which could be active at a time, as: 
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A simple and useful form of penalty functions is given by: 
4 2 
w= ) r h. (ee 
7 
where i stands for violated constraints only 


re is a penalty factor 


h. is the value of constraint violation given by one 


of relations (2.28) - (2.33) 
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Using penalty function (2.34) the objective function is 


given by: 
ME Ea PA (2.35) 


where re is the original objective function. The function f is the 
one to be used in the minimization process. Note that if all 
constraints are satisfied (which is true at the optimal solution), 


the function f is in fact the function ie: 


Although proper minimization requires the solution of a 
sequence of unconstrained minimization problems for monotonically 
increasing res it was observed that for practical power system problems, 
two or three such sequences, each consisting of two or three iterations, 
would be sufficiente**.- Furthermore, as will be shown in this thesis, 
even one such sequence would be enough. However, to do this the 
factor rs Should be chosen in such a way that the increase in the 
value of the objective function due to constraint violations is large 
enough to be sensed by the solution algorithm, but not so large as 
to divert the problem into constraint satisfaction rather than cost 
function minimization. All of this means that in the latter case the 


factor r. should depend on the value of the objective function ee 


2.4 Acceleration 


Although relation (2.21) gives both the direction and the 
size of the minimization step, acceleration would be in order due to 


the fact that the objective function is not really quadratic in system 
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variables, and that the quadratic approximation that Newton's method 
assumes, might not be satisfactory to produce the appropriate step 


size. Relation (2.22) is replaced by: 
=u! +c au’ (2.36) 


where au’ is given by (2.22), and c is an acceleration factor. 


The scheme used here to determine c for any iteration is 
aimed to minimize overshooting and oscillation around the optimal 
solution or a constraint boundary. The factor c is always kept at 
unity until the minimum or a boundary of a functional constraint is 
approached from inside. c is then reduced in such a way as to damp 
any possible oversnoot or oscillation. In more detail, this is 


achieved as follows. 


Before applying the correction Au to the control variables, 
a test is made with c=] to assure that a decrease in the objective 
function would result. If this is not the case, it means that over- 
shooting the solution or constraint violation would result. c is, 
thus, halved and the process is repeated as many times as necessary 
until a decrease in the objective function is achieved or the convergence 


criterion is satisfied. 


To perform the above mentioned test, it would be necessary 
to solve a load flow problem to obtain the value of the x vector 


corresponding to the value u+c Au. Since the value of c is not as 


critical with Newton's method as it would be with other methods, such 
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as steepest descent, it would suffice to solve a linearized load 
flow problem rather than the exact one. The linearized solution 


corresponding to a change au’ in control variables is obtained as: 

pa ce ae) ga (2.37) 
with ax' = S aul (2.38) 
where S is the sensitivity matrix given by: 


] (2.39) 


Note that the first matrix on the right hand side of equation 


(2.39) is the inverse Jacobian which is already available either 
re) 
explicitely or ina factored form. Only the second matrix [ ae 


needs to be computed. Its elements are very simple as shown in 


Appendix A. 


Substitution of (2.39) into (2:38) gives: 

3g 3g 

ie eo er ee 1 
ee (2.40) 
A vector » is then defined as: 


(2.41) 


09 
Note that this product exists since [ uae is a 2N, x (2N oT) rect- 
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angular matrix and au' is a (2N oT) order vector. 


Equation (2.40) can now be written as: 
; 3g 
i = 5 ~| 
i SPL red bh (2.42) 


which is of the form (2.24) with Gen oc (Gada LU arepimesentcuon|) 


a linear repeat solution using the same Jacobian matrix. 


Although the objective of the acceleration scheme is to 
eliminate, as completely as possible, overshooting and oscillations, 
these may still occur due to one or both of the following reasons: 
1) The solution of the linearized load flow problem may not quite 

agree with the exact solution. The former may indicate no over- 
shooting, whereas after applying the latter, at the begining of 
the next iteration, overshooting does occur. Although, this can 
be avoided by using the exact solution throughout, the time 
consumed to obtain one or more such solutions per test would be 
definitely larger than the time consumed in one or two extra 
iterations that may be needed otherwise. Moreover, the disagree- 
ment between the two solutions will be minimal in the final stages 
of the process, due to small Au, when overshooting is more likely 


to happen. 
2) It may happen that overshooting occurs and yet Fey < f.. This 


indicates that the overshoot is small; or, if a constraint violation 


is involved, that r. is small and a solution well outside the 


-_ 7, 7 : ; 


ri Balt i 


} : ; ’ ¥) J mn ’ 


r ap ‘yD 
Te 7 ae ts Lt of) agate hhh 
2 - 7 
| we oy 
- T yi 7 7 
: x b ' - Gc 
oe (> fee 
SA ~~. an \ - a = . 7 KG 
af _ 
g oa « » 
nN - bt Ms = q i: - 
iA PAL 
it oe wet ealnes 1 hii : wa 5% a) “Si ol 


” nko. emy { i) ogee i siiaheee Ki it 


=: 


ine nethratson Sivan WI ieee iy (GLa 497A 
¥ : =f 4 7 rT) 


’ > d i 
: aS (is COLE GMS | ea teeo) sae a ariniog es ote 
’ ; i} ! = 
varia otro WrAd gu oom edlagh "aaa it J ain ont 
i, Ae te il iy q an 4 Ty A ay dq. dal sete “ at wt 


: Gly est hnd vom. secre) Peal pisutae 7008 ha aa 
M : ie) 5 - scan 


; ity iva ni ylagh wade een sontia ode, 
" 
eeits i rin Ag) ie 1h) ni] Th ae si, gt Hibbs 
- = a : 
i lini oo el) ¢ mid neh ee ¥ 
dA Tin jn9 ‘2 i \ Ler 4A ey Th od 4 waren ‘ 
‘ ¥ ; Pe = ) f ' 
eer, ty een , ny SINC wld ALT a Hail WE 
ae 5 ! - 1 a 
op Fetal h. odd cavoeage”, 6iineiongs.? $90 Ale ant Aa ‘Sa y: 
an F 7 t i. ee Leg mt 
7 ' es Tewr?r ata: ty. srl ain od rN wlio tye ao eh woh c 
+) “ee ys ‘fot 


: ¥ Tomi! ia & Ue er NM re ah AEE 2 03 ee " yr " 
| wal ri, oe ie 
nm sy m 7 ay f ) pan 
Cee 

F 


an ae 
.T pa a. vy. ae git 
Manis 
i ens he 
a ob: 
_ 


36 


feasible region will eventually be obtained. 


An optimal acceleration factor (in the sense that it 
minimizes f in the direction Au) may be used resulting in fewer 
iterations. However, the extra time required to solve an associated 
nonlinear load flow problem and the fact that the objective function 
should, for the purpose of computing this optimal acceleration factor, 
be approximated to a second or third degree polynomial in c, may 
lead to almost the same convergence behaviour. Add to that, that the 
second reason above is general and applies to any acceleration 
mechanism. Therefore, the extra complexities in programming such 


optimal acceleration scheme do not appear to be justified. 


2.5 Voltage Controlled Busses 


A voltage controlled bus is a bus at which the voltage 
magnitude is controlled through the introduction of reactive power 
generation. As a result, voltage magnitude should be considered as 
a specified or control variable, and reactive power at the bus is 
left to change within the operating limits of the generation equipment. 
This means that a number of variables (voltage magnitude) corresponding 
to the number of such busses will be moved from the dependent variable 
vector x to the control variable vector u. Also introduced at each 
such bus is a functional constraint of the form (2.14) defining the 
range in which reactive generation at that node can change. Since the 
voltage magnitude is no longer a dependent variable and reactive 
power ceases to be a fixed amount, the reactive power equation 


corresponding to each of these nodes is removed from load flow equations 
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In other words, a voltage controlled bus is to be treated 
as a generator (or swing) node as far as reactive power is concerned, 


and as a load node as far as real power is concerned. 


The order of the different vectors in cases involving No 
such nodes will be abel + N for the control variables u, and 
EN ONG for each of the vector g representing the load flow equations 


and the vector x of the dependent variables. 


2.6 Qualitative Evaluation 


The aim of the solution method is to reduce computer 
storage requirement, and to achieve better convergence characteristics 
than other methods. Based on these two objectives, a qualitative 
evaluation of the method can be carried out. The methods of 
Dommel and Parana evel and Sasson et 31 l22 1] were considered for 
comparison purposes. Also given is a rough estimation of computation 


time requirements as compared to these two methods. 


2.6.1 Computer Storage 


It is customary, when one speaks of computer storage, to 
consider large systems, and estimate the amount of storage required 
for non-zero elements only. Since different matrices differ in their 
sparsity characteristics, and even the sparsity of the same matrix 
differs from one system to another, it would be appropriate when 
comparing storage requirements of different methods to use the concept 


of a full matrix and then to comment on what effects sparsity will have 


on each. 
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Apart from storage requirement of the vectors x, u and 
other common quantities as power generations and loads, which should 
be the same for all methods, computer storage requirements of the 
developed method and the other two are shown in Figure 2.4 and Table 
2.1. For simplicity no voltage controlled busses are assumed to be 
present. Table 2.2 also gives the amount of storage required for 
five standard test systems. Voltage controlled busses were considered 


in thus case. 


As can be seen, the method of Sasson et 41 l22] requires 
the most storage for their Hessian matrix. Dommel and Tinney's 
method may require more or less storage depending on the size of the 
system and the number of its generator busses. Most of their 
storage area is for the Jacobian matrix, however. The proposed method 
requires the least storage of the three methods. Moreover, as in 
the case of Dommel and Tinney's method, most of the required storage 


is for the Jacobian matrix. 


As for sparsity, it has been shown that for typical power 
systems the Hessian matrices are more full than the Jacobian matrices*. 
Therefore, sparsity techniques of storage will give more advantage to 
those methods where the Jacobian occupies most of the storage. Hence, 
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storage savings effected by the proposed method can be emphasized even 


further. 


2.6.2 Convergence 


Since the iteration cycle of the proposed method and that 
of Dommel and Tinney are identical, convergence behaviour as 
judged by the number of iterations required to obtain the solution, 
will depend on the way the control variables are modified. Instantly, 
second order gradient methods, as used in this thesis, emerge as the 
Superior compared to the steepest descent correction as used in 
Dommel and Tinney's algorithm. Moreover, that latter method is so 
sensitive to the acceleration factor c, that it will not converge in 
a reasonable number of iterations, or at all, unless the factor c is 
carefully chosen. This type of sensitivity does not exist in Newton's 


method. 


Sasson et al used Newton's formula as used here. One then 
should expect that the two methods will have the same convergence 
characteristics. However, Since in their method a search of the whole 
space rather than the feasible region is required, more iterations 


will be needed to arrive at the solution. 


2.6.3 Computation Time 

Accurate assessment of the proposed and the two other methods as 
far as the computation time needed to obtain a solution is possible 
only if programs of compatible efficiency were written for each 


algorithm and run on the same computer. However, a rough comparison 
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is still possible in the following way: 


1) While the proposed method uses the Hessian matrix of the second 
order partial derivatives of the objective function with respect 
to the control variables to compute the corrections in these 
control variables, Dommel and Tinney's method+ |» uses only the 
gradient. However, as mentioned before, this Hessian matrix is 
of low order, extremely sparse and its elements are very easy to 
compute. Moreover, the load flow problem, which is to be solved 
in both methods, is of much lower order in the proposed method than 
in Dommel and Tinney's, thus requiring less computer time. 
Although part of this time saving will definitely be used in 
handling the Hessian, in no case will total time per iteration 
required by the proposed method exceed that required by Dommel 


and Tinney's method. 


2) Although the proposed method and that of Dommel and Tinney require 
load flow solution each iteration (such solutions usually require 
one or two Newton iterations), Sasson et 4122] uses a much larger 
and less sparse matrix, and, thus, will not offer any savings in 


the time required per iteration. 


One can now estimate that the three methods will require 
about the same amount of time per iteration. The total time required 
to obtain the solution will then depend on the number of iterations 
needed to produce that solution. Here the decision will be in favour 


of the proposed.method. 
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CHAPTER III 
OPTIMAL SOLUTIONS OF SMALL SYSTEMS 


This chapter deals with the application of the method 
developed in the last chapter to small power systems. In handling 
these small systems one does not have to worry about some of the 
problems encountered with large systems. For instance, the matrices 
involved are of low order and proper matrix inversion, for which ready 
routines are available, rather than factorization, can be used. 
Furthermore, storage requirements are very limited resulting in that 
the natural order of nodes would be satisfactory, and an optimal 
ordering scheme can be omitted. The result is a Simple computer 


code and quick assessment of the success or failure of the method. 


[32] 


Two standard test power systems were studied, a 5-bus system 
and the IEEE 14-bus test system. Data of the two systems is given in 
Appendix B. Four optimization problems were solved for each system 
using the University of Alberta IBM 360/67 Computer. The problems are 
minimum operating cost, minimum losses, minimum fuel comsumption, and 
combined cost-fuel minimization. The two latter problems are defined 


and formulated in this thesis for the first time. 


3.1 Minimum Operating Cost Problem 
This is the exact economic dispatch problem. The objective 
function is the total cost of generation, which is assumed to be of 


the form: 
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The constants ass b. and C. for both test systems are given 
in Table 3.1. Since, in the method, generator and load nodes are 
handled separately, the natural order of nodes in both systems is 
changed so that generator nodes are numbered first, voltage controlled 


busses follow, and at the end come load busses. 


Tables 3.2 and 3.3 give the solution of the problem , for 
the 5-bus and 14-bus systems respectively, when voltage magnitudes 
at generator and voltage controlled busses were kept fixed. In the 
case of the 5-bus system a flat start was used. However, in the case 
of the 14-bus system, a different starting point was chosen, since 
the flat start proved to be too far from the optimal solution, resulting 
in a slow convergence. This behaviour is typical of any iterative 
method. It did show, however, that the algorithm can easily handle 
violated functional constraints. The operating point, which was outside 
the feasible region for the flat start, was brought inside in two 


iterations. After that convergence was slow. 


Since the optimal solutions were inside the feasible region, 
and to test the algorithm's ability to handle functional constraints 
if they are violated near the solution, which is a more sensitive 
problem, several runs were carried out on the 5-bus system. The 


results of these runs are given in Tables 3.4 - 3.8. 


In Table 3.4, the voltages at generator busses no. 1] and 2 
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table 3.1 Cost Coefficients. of 5- and 14- Bus 


Systems 


Total Cost of Generator = 760.95 $/Hr. 


Table 3.2 Optimal Solution of Minimum Cost Problems for the 5-Bus 


System (fixed Generator Voltages) 
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Table 3.3 Optimal Solution of Minimum Cost Problem 
for the 14-Bus System 


(Fixed Generator Voltages) 
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f = Total Cost of Generation = 757.78 $/H. 


* Constraint violation = 9x107° peur 


Table 3.4 Optimal Solution of Minimum Cost Problem for the 5- 


Bus System (Free Generator Voltages) 
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were free to change within the limits shown. In this case the 
unconstrained optimal solution was outside the feasible region resulting 
in a voltage constraint violation at load node no. 5. However, the 
algorithm was able to locate the constrained solution holding the 
voltage magnitude at that node within 9x10" p.u. Of the constraint 
boundary. As a result generator voltages could not rise to their 


maximum limit as one expects in such a case. 


In Tables 3.5 - 3.8, node no. 3 was changed into a voltage 
controlled bus ener reactive generation was introduced and the 
voltage was allowed to change within the limits shown. The intention 
is to show how nonlinear functional constraints, such as real or 
reactive power limits, are handled if violated. Two ranges of reactive 


generation at node no. 3 were assumed: -0.5 to + 0.5, and - 0.4 to + 0.4. 


Tables 3.5 and 3.6 give, respectively, the solutions of 
these two cases when voltage magnitudes of generator nodes were kept 
fixed. In the first case, the solution is well inside the feasible 
region, and it represents the unconstrained (as far as functional 
constraints are concerned) minimum of the objective function. Although 
this unconstrained solution lies outside the feasible region of the 
other case, a solution was obtained at which reactive generation at 
node no. 3 exceeded its maximum limit by only > 5ExI0N, which shows 


the effective handling of this type of constraint. 


Tables 3.7 and 3.8 give the results for the same cases when 


voltage magnitudes at generator nodes no. 1 and 2 were allowed to 


change as well. Although in the first case, reactive power generation 
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f ="Total Cost of Generation = 757.57 S/Hr. 


Table 3.5 Optimal Solution of Minimum Cost Problem for Modified 


5-Bus System (Fixed Generator Voltages) 
-0.9<Q.<0.5 


f = Total Cost of Generation = 757.65 $/Hr. 


: -3 
* Constraint violation = 2.55xl0 ~ p.u. 


Table 3.6 Optimal Solution of Minimum Cost Problem for Modified 


5-Bus System (Fixed Generator Voltages) 


-0.4<0<0.4 


J 


f = Total Cost of Generation = 755. 


* No Constraint violation 

Table 3.8 Optimal Solution of Minimum Cost Problem for Modified 
5-Bus System (Free Generator Voltages) 
-0.4 < Q, < 0.4 
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at node no. 3 stayed well within its limits, the voltage magnitude 


oh In 


at node no. 5 exceeded its maximum limit by only 5.19x10~ 
the second case, in which the constraints on reactive generation at 
node no. 3 and voltage magnitude at node no. 5 were violated during 
the minimization process, the solution obtained was just inside these 


constraint boundaries. 
3.2 Minimum Loss Problem 


This is actually the problem of minimizing total real power 
generation nike TE determined by the total load plus system losses. 
The problem is often termed "optimal reactive power flow" because 
it is mainly the reactive flow in the transmission system that 
determines the losses. For lower losses, reactive generations should 


be close to the loads. 


The objective function is given by: 
tn) i, (322) 


The solutions for the 5- and 14-bus systems are given in 
Tables 3.9 and 3.10 respectively. In the 5-bus system the only 
constraint violation was in the voltage magnitude at bus no. 5 which 
exceeded its maximum limit by only 7.85x10"" p.u. thus holding 
generator voltages from reaching their maximum limits where the 


unconstrained (with respect to functional constraints) minimum is 


located. In the 14-bus system, real generations at generators 1 and 2 
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f = Total System Losses = 0.0438 p.u. 


* Constraint violation = 7.85x10~" p.u. 


Table 3.9 Solution of Minimum Loss Problem for the 5-Bus 


System (Free Generator Voltages) 
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f = Total System Losses = 0.06678 p.u. 


* Constraint violation = eGo Cie pau, &* Cons trainteviolation SB 10x10) pel 


Table 3.10 Solution of Minimum Loss Problems for 
the 14- Bus System 


(Fixed Generator Voltages) 
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exceeded their limits. However, these constraint violations were 
successfully controlled during the minimization process, such that 


at the solution, maximum violation was 7a56R) OG pi. Ute 
3.3 Minimum Fuel Problem 


Although generation cost at a power plant consists mostly 
of fuel cost, it includes the cost of labour, maintenance, etc. 
This portion is not the same for all plants within a system. For 
instance, an old plant needs maintenance more frequently than a new 
one; a gas turbine power plant needs a different amount of labour 
than a steam plant of the same capacity. Furthermore, fuel transportation 
costs are different for different plants. Therefore, generation 
cost functions are not indicative of the true fuel consumption in 
a given system, and scheduling generation according to a criterion of 
minimum generation cost does not necessarily mean that fuel consumption 


is also at a minimum. 


The importance of the minimum fuel problem evolves from 
environmental and resource conservation grounds. On one hand oxides 
emission will be held as low as possible , and on the other hand, by 
consuming a minimum amount of fuel for a given load condition, it is 
possible to prolong the life time of depleting fuel resources. The 
: In fact oxidesemission does not depend only on the amount of fuel 


burned but also on the type of fuel and how perfectly it is burned. 
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idea is to make use of efficient power plants to produce most of 
the load requirement irrespective of the cost of such production. 
Thus, in return for the above mentioned advantages one is to pay a 


premium in the form of more expensive energy. 


Let oF be the proportion of non fuel costs in the generation 
cost function of generator i. Also assume that generator i receives 
LESMiUelodteda COS t.0f Y; dollars/million B.T.U.'s. If generation cost 
is assumed to be of the form (3.1), then, fuel consumption of 


generator 7 in million B.T.U./Hr. is given by 


l-a 5 
f. = Y; (a0 b. P. + c, P. ) (253) 
The objective function in this case will be 
faz ptt 
j 
I-a, 2 
= d (a, + b. P. + C; Phe) (3.4) 
=e (amie gbvbaP we-rnceneP aS) (325) 
1 j 7] 


The problem was solved for the 5- and 14-bus systems. The 
coefficients oF and y; are assumed as given in Table 3.11. The 
coefficients as» b. and C. are kept as given in Table 3.1. Also shown 
in Table 3.11 are the coefficients a,'s b.' and Cio of the objective 
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Table 3.11 Coefficients of the Objective Function of 
the Minimum Fuel Problem for the 5- and 


14-Bus Systems 


Tables 3.12 and 3.13 give the solution for the 5-bus 
system for fixed and free generator voltages respectively. Table 3.14 


gives the solution for the 14-bus system. 
3.4 Combined Cost-Fuel Minimization 


Each of the minimum cost and minimum fuel problems is based 
on one global minimization criterion governing the whole system, 
consequently the solution will satisfy this criterion globally rather 
than at each individual generating plant. Sometimes, that is not 
quite what is desired. For instance, in addition to scheduling 
generation to achieve global minimum cost, one plant, a group of plants, 
or all plants may be required to consume a minimum amount of fuel to 
overcome a fuel shortage condition, or to achieve a reduction in 
pollution level in some locality. This means that scheduling of 
generation would be based on two criteria. In other words, a compromise 
should be arrived at which may not quite satisfy each criterion 


individually, but yet it minimizes their sum. One of these two criteria 
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f = Total Fuel Consumption = 1319.00 Million BUT/Hr. 


Table 3.12 Solution of Minimum. Fuel Problem for 5-Bus System 


(Fixed Generator Voltages) 


f = Total Fuel Consumption = 1314.27 Million BTU/Hr. 


* Constraint violation = 8.12 x 107 peur 


Table 3.13 Solution of Minimum Fuel Problem for 


5-Bus System (Free Generator Voltages) 
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f = Total Fuel Consumption = 2117.86 Million BTU/hr. 


Table 3.14 Solution of Minimum Fuel Problem for 14- 


Bus System (Fixed Generator Voltages ) 
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should be global while the other may be only of local concern. 


The objective function at a generating node would, then, 


be: 


foes O. FB oats (376) 


where ie and uF are the cost and fuel functions at node i, and 
c P 
oF and Oo; are priority factors. 
Cc fi 


It is important that both functions be in the same units 
otherwise one will dominate the other (see Tables 3.1 and 3.10). One 
possibility is to transform generation cost into fuel units by dividing 
the generation cost function at each generator by fuel price at that 
node. However, since cost minimization is usually of global rather 
than local interest, and, it is more likely that fuel minimization is 
of local concern only, it is more appropriate to convert fuel 
consumption into dollars/Hr. Merely multiplying fuel functions by 
corresponding fuel prices will not be satisfactory. The resulting 
functions in this case will not be indicative of relative fuel 
consumption at system plants, since such prices depend on the type of 
fuel and include different transportation and handling costs. Thus, a 


base price should be used instead. The functions f, will, then, be: 
iP 


fa 
(a, +e DH IPAME ECE P. ) 1; (327) 
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where Yp is the selected base fuel price. 


The total objective function in this case will, then, be 


given by: 


(3.8) 


Diseatea plant, Faecertainecri terion) isso funor; concerns its 


priority factor is set to zero. 


Based on a base fuel price of 0.40 $/million B.T.U. cost- 
fuel function coefficients for the 5- and 14-bus systems are given in 
Table 3.15. It should be stressed here that these coefficients are 
used only when both criteria are active at a given node. Otherwise, 


the coefficients of Table 3.1 should be used. 


Table 3.15 Cost-Fuel Coefficients of the 5- 


and 14-Bus Systems 
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Several runs were carried out for different situations. 

It is assumed that whenever a criterion is active at a plant, it is 

considered to be as impertant as the global criterion, i.e. priority 
factors are taken as unity when their corresponding criteria are of 

concern. 

Two situations were considered for the 5-bus system. Two 
cases corresponding to fixed and free generator voltages were run 
for each. In the first situation, it is assumed that the objective 
is to achieve an overall economy in both generation cost and fuel 
consumption. The coefficients of the objective function would, then, 
be those of Table 3.15. The solutions for this situation are given 


TnmbdvteseoelOrand o.1/. 


In the second situation, what is required is to achieve an 
overall economy in generation cost such that economy in fuel consumption 
of generator no. 1 is also realized. Therefore, objective function 
coefficients for generator no. 1 would be drawn from Table 3.15, while 
those for generator no. 2 from Table 3.1. The solutions for this case 
are given in Tables 3.18 and 3.19. Observe how generator no. 2 has 
picked its maximum share of the load to allow the minimization of fuel 


consumption of generator no. 1. 


In the case of the 14-bus system, three different objectives 
were considered. The first is the same as the first case of the 5- 
bus system, i.e. an overall economy in cost and fuel is desired. Table 


3.20 gives its solution. The second is to achieve an overall cost 
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(Cost and Fuel Consumption) = 1292.68 $/Hr 
Cost of Generation 702.7 Oe oy te 
Fuel Consumption 1324.78 Mitlion, BIU/HY: 


Table 3.16 Solution of Minimum Cost-Fuel Problem for the 


5-Bus System (Fixed Generator Voltages) 


0.376 0.6 | 
0.400] 0.6 
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Total (Cost and Fuel Consumption) 
Total Cost of Generation 
Total Fuel Consumption 
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* Constraint violation = 1.37 x 107° p.U. 


Table 3.17 Solution of Minimum Cost-Fuel Problem for the 


5-Bus System (Free Generator Voltages) 
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56.76 $/Hr. 
94.58 $?hr. 
3 Million BTU/Hr. 


Total Cost and Fuel Consumption at Gen. no. 1. 
Total Cost of Generation 


= 9 
= 7 
tie = Fuel Consumption of Gen. no. 1 = 405.4 
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Table 3.18 Solution of Minimum Total Cost and Fuel Consumption 
of Generator no. 1 for the 5-Bus System 


(Fixed Generator Voltages) 


f = Total Cost and Fuel Consumption of Gen. no. 1 = 951.36 $/Hr 

Te = Total Cost of Generation = Pe CAO sy lete 

Teel = Fuel Consumption of Gen. no. 1 = 395.37 Million BTU/Hr 
| 3 


* no constraint violation ** constraint violation = 5.77 x 10 ~ p.u. 
Table 3.19 Solution of Minimum Total Cost and Fuel Consumption 
of Generator no. 1 for the 5-Bus System 


(Free Generator Voltages) 
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1991.42 $/Hr. 


Total (Cost and Fuel Consumption ) = 
Total Cost of Generation = Cove fey bine 


Total Fuel Consumption 2128.88 Million BTU/Hr 


Table 3.20 Solution of Minimum Cost-Fuel Problem for 
the 14-Bus System. 


(Fixed Generator Voltages) 
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Total Cost and Fuel Consumption of Gen. 1 &2 = 1709.28 $/Hr. 
Total Cost of Generation Wile CCRS HY. 


Fuel Consumption of Generators 1 & 2 1345765 Mild ion BTU/hr. 


* Constraint violation = 1.099x10" Deu. 

Table 3.21 Solution of Minimum Total Cost and Fuel 
Consumption of Generators no. 1 and 2 for the 
14-Bus System 


(Fixed Generator Voltages) 
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Cost and Fuel Consumption of Gen. 2 = 1241.19 $/Hr. 


Total Cost of Generation 1157.59 $/Hr. 


Fuel Consumption of Generator 2 209.01 Million BTU/Hr. 


* Constraint violation = 1.385 x ios D.uU. 


Table 3.22 Solution of Minimum Cost and Fuel Consumption 
of Generator no. 2 for the 14-Bus System 


(Fixed Generator Voltages) 
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economy and also to economize in fuel consumption of generators 

no. | and 2. The intention was to show that generator no. 3, whose 

fuel consumption is of no interest, will provide its maximum generation, 
and the rest of the load will be scheduled between generators 1 and 2. 
ihe solution*of this *case=is given’ in Table 3.21le the last case 1s 
that overall cost economy is sought such that at the same time fuel 
consumption of only generator 2 is economized. The purpose here is to 
show that generator 2 will provide its minimum generation and the 

rest of the load will be scheduled by generators 1 and 3, where fuel 


consumption is of no concern. Table 3.22 gives the solution. 


3.5 Convergence Behaviour 


Tables 3.23 and 3.24 show the number of iterations needed 
to obtain the solution of the cases studied and reported in the previous 
sections, for the 5- and 14-bus systems, respectively. Also shown is 
the maximum violation in functional inequality constraints. Parameter 
constraints were strictly observed by virtue of equation (2.26). 
Equality constraints were considered to be satisfied when power mismatch 


of iar’ or less was obtained in the load flow portion of the problem. 


In all cases of the 5-bus system, a flat starting point was 
used. However, in all cases of the 14-bus system it was too far from 
the optimal operating conditions resulting in slow convergence (e.g. 32 
iterations in the minimum cost problem). However, this does not detract 
from the effectiveness of the method. The flat start is not a realistic 


operating point as all generators are in phase with the reference 
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Table 3.23 Convergence Characteristics for the 5-Bus 
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Sasson et. al 


Proposed Method 
ter. Max.Viol. 


(Cost and Fuel) 


Min. (Cost and Fuel of Gen. ike 099x107° 
Tah 2) 
= 1385x1072 


Table 3.24 Convergence Characteristics for the 14- 
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generator. It also does not take into consideration the available 
information and the experience gained on the system under study, and 
it is likely that a different starting point will be chosen based on 
practical judgement of what value each variable is likely to assume. 
Even with the lack of this information the results of preliminary 


studies can be used. 


Recognizing that, it can be seen that the method possesses 
an excellent convergence rate. Furthermore, this rate does not depend 
On the size of the system, since for both systems, solutions are 


obtained in a comparable number of iterations. 


It should be mentioned here that in the cases which required 
more than 6 or 7 iterations, the objective function was almost flat 
near the solution, and a large number of iterations were used for 
little improvement. The minimization process would have terminated 


much earlier if the improvement was taken as the convergence test. 


Tables 3.23 and 3.24 also show the number of iterations 
required to obtain the solution of the minimum cost problem, and their 
associated maximum constraint violations, by the method of Sasson et al 


This is the only problem solved in their paper. Since in both cases, 


the solution is well inside the feasible region, the violations are 
in the equality constraints. These figures are practically unaccept- 


able and a few extra iterations are needed to improve on them. 


Although exact comparison between the two methods is not 


possible with the lack of the same starting point*, those figures 


a 


* no starting point is mentioned in ref. [22]. 
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are a rough indication of the relative effectiveness of the methods. 
3.6 Solution Qualification 


By the term "solution qualification" is meant the invest- 
igation of whether or not a solution of a problem does meet the 
problem requirements. To achieve this qualification, one should, first 
investigate the power schedule of each individual solution to see if 
it agrees with simple intuition provided by the problem itself. 
Secondly, a comparative study of the solutions of different problems 
in terms of the value of the different objective functions is carried out. 
Finally, one is to show that the solutions obtained are indeed the optimum 
in their respective cases. It should be mentioned that, from the 
point of view of constraint satisfaction, all solutions, despite the 
violations involved in some, are quite acceptable for all practical 
purposes. Consequently, each solution is considered to satisfy its 


constraint requirements. 


Tables 3.25 and 3.26 give power generation schedule, and 
the values of the different objective functions involved in each case 
for the 5- and the 14- bus systems, respectively. For easy reference 
the different cases were numbered, though their numbers do not agree 
with the order of their presentation in this chapter. Also included 
in Table 3.25(a) is the solution of the minimum loss problems of the 
5- bus system for fixed generator voltages, which were extremely close 


to the flat starting point. 


Of all cases studied on the two test systems, power generation 
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ait} 


0.05167 
0.05021 
0.05083 
0.05756 


0.05008 
0.04471 
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Losses 
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0.04381 
0.04469 
0.05171 


0.04379 
0.03874 


0.03868 


Fue] 
106 
BTU/Hr. 


1335299 
1324.68 
1319.00 
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1329259 
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106 
BTU/Hr. 
1334730 
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1314.27 
1327228 


1318773 
1328.42 
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Fuel of 
Gen. | 


Table 3.25 Power Schedule and Value of Objective 
Functions of the 5-Bus System 
(a) Fixed Generator Voltages 


(b) Free Generator Voltages 
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schedule of only a few can be readily predicted from the statement 
of the problem. In all other cases, however, a study of the objective 
function would produce an idea of how generation is shared by the 


various generators of the system. 


The first category consists of all cases where local 
criteria were invoked. These are cases no. 4 and 11 of the 5- bus 
system and cases 4 and 5 of the 14- bus system (see Tables 3.25 and 
3.26). In such cases one expects that a particular generator (or 
a group of generators) will provide minimum generation when it is 
(they are) required to consume a minimum amount of fuel, provided 
that the balance of the load can be absorbed by the other generators 
of the system. If this is not the case, those other generators will 
provide their maximum share, and the constrained generator(s) must 
absorb the balance. Investigation of Tables 3.25 and 3.26 shows that 


this is exactly what has been achieved in the above mentioned cases. 


The second category comprises all other cases where minimi- 
zation criteria were global. Although generation schedules in these 
cases are less predictable, one should be able to predict the way 
generation will be shared between various generators by studying their 
respective contribution to the objective function. A generator 
contributing less to the objective function will definitely carry 
a largershare of the load. This is exactly what happened in all cases 
of this category. Even in case no. 5 of the 14- bus system, generators 
1 and 3 shared the balance of the load in such a way as to minimize 


the total cost of generation, the global criterion. 
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Investigation of the values of the total cost of generation, 
system losses, and total fuel consumption as given in Tables 3.25 
and 3.26 shows that any of these quantities are minimized at the expense 
of the two others. A compromise solution can be obtained if all 
participate in the objective function, as was in the case of the cost- 


fuel problem. 


Invoking the local minimization criteria resulted in a 
considerable increase in the cost of generation for both systems, and 
also in fuel consumption in one of the 14-bus system cases. The reason 
is that generation was forced to shift to the more expensive source, 
so that the local criterion can be observed. It should be mentioned 
here that the choice of these local criteria was intended to show such 
effect. This may or may not be the case in a realistic situation 
depending on the cost functions associated with the various generators 


of the system. 


Investigation of Table 3.25 will also show that the solutions 
did not fail in showing the well known effects of increasing the 
voltage level in a system, and those of introducing reactive generation. 
Observe that generation cost, system losses and fuel consumption are 
lower in the cases of free generator voltages than those of fixed 
generator voltage. Also, they are lower in the cases where there is 
reactive generation at bus no. 3 than the cases which do not have such 


generation. 


What is left now is to show that these solutions are indeed 
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the optimal solutions and not just a boundary point. First, invest- 
igating Tables 3.25 and 3.26 shows that the minimum of each objective 
function, local or global, does occur opposite the case requiring 


the minimization of such an objective. 


Secondly, the convergence behaviour, as indicated by the 
value of the acceleration factor, points out that these solutions are 
the optimal. The acceleration factor in most cases stayed at unity, 
while in the rest it decreased gradually when the solution was 
approached. Should it change sharply it would be an indication that 
a boundary point, rather than the optimal, will be obtained. This will 


be discussed later. 


Thirdly, Figure 3.1 shows the equality constraint curve of 
the 5- bus system, with fixed generator voltages, on the P1-Po plane. 
Also shown are some constant total generation cost contours. The 
diagram represents the minimum generation cost problem. All operating 
points of the system, including those resulting during the minimization 
process, must lie on the equality constraint curve. This curve will 
cross many constant cost contours representing different cost. However, 
it is the point at which the curve is tangent to one of these contours 
that represents the optimal solution. The optimal operating cost will 
be that represented by that contour. yes was indeed the case. For 
instance the contour corresponding to a cost of 700 $/Hr. represents 
a lower cost than the optimal. However, any point on that curve does 


not satisfy the equality constraints, thus, that cost can not be realized. 
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Figure 3.1 Equality Constraints and Constant Cost Contours 
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The lowest realizable cost is that represented by the 761 $/Hr. contour 
which touches the curve of the equality constraints at the optimal 
solution. Higher cost contours will each intersect with that curve 


at two points resulting in suboptimal solutions. 


For a lossy system the curve of equality constraints will 
be concave upwards ,whereas the constant cost contours are concave 
downwards. This will ensure that a single minimum occurs. In this 
case, however, since system losses are small (typical for power 
systems) and there is little difference between the light load costs 
at the two generators, this curvature is small. The result is that 
the true minimum point is difficult to locate with any realistic 
convergence criterion. However, in practice, this is unimportant as 
the cost (or the operating point depending on the convergence criterion) 


will be within the preset tolerance of the actual minimum. 


Similar diagrams can be drawn for the other objective 
functions of this case*. However, although the cases of the 5-bus 


system with free generator voltages, and those of the 14- bus system 


ac A 


* In the minimum loss problem, the constant loss contours will be 
discontinuous, and each will reduce to two points bracketing the 
optimal point. This is because system losses are defined only 


on the curve of equality constraints. 
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follow the same theory, it is difficult to draw such diagrams as they 


are multidimensional. 


Finally, other methodst2! »22,33] have treated two of the 
problems presented in this chapter, namely, the minimum cost problem 
of the 5- and 14- bus systems. For all practical purposes their 


solutions agree with the solutions obtained by the proposed method. 


8] 


CHAPTER IV 
OPTIMAL ORDERING OF LARGE SYSTEMS 


As with any solution method, in applying the proposed 
method to a large system, one has to deal with several problems 
that are not involved in the case of small systems. For instance, the 
computer storage required to store the system matrix would be 
prohibitive and the efforts of inverting such a matrix very expensive. 
However, as mentioned before, these problems have been greatly eased 
in the case of power systems due to the extreme sparsity of the 
matrices involved, e.g. the Jacobian matrix. Only non zero elements 
of the matrix would be stored, but then, one has to keep track of 
their positions in the original matrix. Matrix inversion as such 
must be avoided leading to the use of an elimination process (e.g. 
Gaussian elimination) to triangulize the matrix and then using back 
substitution to obtain the solution. An optimal ordering scheme 
to order the rows of such a matrix would be in order to preserve its 
sparsity during the elimination process, thus reducing the extra 


storage that would be needed. 


In this chapter a new method of optimal ordering is 
described, and compared to two other established methods. Three 
standard test systems are used. They are the IEEE 30-, 5/-, and 
118- bus test systems. Data of these systems is given in Appendix 
B. Optimal cost solution of the 30-bus system is also presented as a 


representative of optimal solutions of large systems. 
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4.1 Optimal Ordering 


An absolute optimal ordering scheme of the elimination 
process, through reordering system nodes, would result in the 
least possible terms in the table of factors*. An efficient 
algorithm for determining such absolute optimal order has not been 
developed, and it appears to be a practical impossibility. However, 
several effective schemes to determine a near-optimal order have 
been developed. Although these schemes do not produce the true 
optimal order, they are known, and hereafter are referred to, as 
optimal ordering schemes. An extensive comparative study of these 


schemes is found in reference 34. 


There are two basic approaches to optimal ordering. The 
first is preordering, and it amounts to renumbering system nodes, i.e. 
the rows and columns of the system matrix, according to the order 
required in the elimination process. Solutions can then be obtained 
by choosing successive pivots down the leading diagonal. The main 
advantages of such schemes are that they are simple to program 
(consequently they are fast to execute),and that the only information 


needed to perform them is a node-branch connection list. 


The best known ordering scheme using this approach is a one 
in which system nodes are numbered starting with that having the 
* The table of factors refers to the triangulized matrix when 


forward operations are stored in the lower triangle. 
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fewest number of connected branches, and ending with that having 
the most. It does not, however, take into account anything that 


[26,27] This scheme will 


happens during the elimination process 
be referred to as scheme I for the purpose of comparison later in 


the chapter. 


The second approach is known as dynamic ordering. System 
nodes are ordered according to the way the table of factors develops 
during the elimination process, rather than from the structural 
properties of the original system matrix. Programming of such 
methods is, therefore, more complicated, as ordering takes place 
during the actual elimination process or at least using its 
simulation. Execution of such schemes is slower than preordering 
schemes; however, they do produce an order closer to the absolute 


optimal. 


Of the dynamic ordering schemes, two are the most familiar. 
In the first, referred to as scheme II, at each step, the next node 
to be eliminated is the one with the least number of connected lines 
in the reduced network i.e. the next variable to be eliminated is 
that associated with the row and column of the reduced matrix 
containing the fewest non-zero off diagonal elements. In the second 
scheme, the next node to be eliminated is the one which will 
introduce the fewest new equivalent lines. This scheme not only 
requires the simulation of the elimination process according to the 


natural order of the remaining nodes, but of every feasible 
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alternative at each step. This scheme is referred to as scheme 


ED. 


The choice of scheme is a trade-off between speed of 
execution and the number of times the result is to be used. Scheme 
I, for instance, is good for problems requiring a single solution 
with no iterations. Scheme III, on the other hand, is good for 
problems requiring a large number of iterations to justify the 
extra time used in its execution. For load flow studies, scheme 


II was found to be the best of the threeL26>27 230534] 


4.2 The Proposed Method 


The main disadvantage of preordering schemes is that 
they do not take into account the changes that occur in the table 
of factors during the elimination process. This makes such schemes 
unsuitable for iterative solutions. On the other hand, dynamic 
ordering schemes require longer time to execute which makes them 


relatively unsuited to solutions that require only few iterations. 


The proposed method falls within the category of pre- 
ordering schemes. However, ordering is carried out in such a way 
to reduce, as much as possible, the structural changes in the table 


of factors, thus reducing the number of the new elements introduced 
during the elimination process. 

It would be instructive to show that in a Sparse symmetric, 
or incident symmetric, matrix, if the rows are ordered such that 


off-diagonal elements become more and more dense aS One moves down 
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along the major diagonal, that the chance of introducing new 

elements during the elimination process becomes more and more 
unlikely. A simple reason is that very few Openings are left for 
fill-in terms. The elimination process will merely modify already 
existing elements. In power systems, each row has at least one 
off-diagonal element*, and it is unavoidable for such terms to appear 
in every row. The strategy will, then, be to restrict these 

elements to the triangular portion below the minor diagonal, and, 

as before, to make them more concentrated in the lower right hand 
corner of the matrix. This would be the ideal structure for the 


elimination process in such casesee/ 23°). 


An ordering scheme of system nodes to achieve such a 
structure is given as follows: 

1) If there are any nodeswith no connected branches, they are 
numbered first. As a matter of fact these nodes will not affect 
the elimination process, and can take any order. However, they 
are numbered first for simplicity. 

2) Numbered next is the node connected to the least number of 
branches. If more than one qualify, the choice could be either 
arbitrary, or according to a descending order of the number of 
branches connected to those other end nodes. 

* For the present formulation of the optimal load-flow problem, it could 

happen that a node, particularly a generator node, may not have any 

connected branch leading to another node of the same set of nodes, 


because load nodes and generator nodes are treated separately. 
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3) The node (or nodes) which is (are) connected to the node just 
ordered, and which is (are) not previously numbered, is (are), 
then given the highest order(s) available, according to a 
descending order of the number of branches connected to each. 

4) Of the remaining nodes, the node having the least number of 
branches is picked up. If no more than one qualify, steps 2) 
and 3) are repeated. However, if more than one are found, they 
are searched for a node connected to the highest ordered node. 
If one is found, it is given the lowest available order, and 
step 3) is repeated. If no one is found, the search is repeated 
for a connection with the second highest ordered node, and so 
on. When these highest ordered nodes are covered, steps 2) 
and 3) are repeated. The process is continued until all nodes 


have been numbered. 


It can be seen that once a low order node is numbered, 
all its connecting nodes are moved away to the minor diagonal 
and beyond. Furthermore, such connecting nodes are numbered such 
that the node having the highest number of connected branches is 
given the highest order. This will ensure that the larger the 
number of the off-diagonal elements, the lower they will appear 


in the matrix being ordered. 


For relatively small systems, the method may not work as 
well as described. The number of nodes which are connected to only 
one other node is quite small. Some of these nodes may be connected 


to the same node and the chances are that one may have to order a 
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node connected to two ore more nodes of which nonehas been previously 
numbered. The result will be quite a few off-diagonal elements 
appearing to the left of the minor diagonal. The Opposite is true 
for large systems. The number of those nodeswhich are connected 

to only one node is relatively large and scattered around the system, 
such that when nodes connected to two or more nodes have to be 
ordered some of these connecting nodes would have been already 
numbered and the off-diagonal elements will be located about the 
minor diagonal and to its right. The number of the off-diagonal 
elements to the left of the minor diagonal will be very small, 


however. 


OUST SaaiuleliUSutea teins Maula GU, Cmca yemmelilnie ave area 
represents the area of the ordered matrix where off-diagonal elements 
will appear. For large systems, an area to the right of the minor 
diagonal will have no elements. It corresponds to those nodes with 
only one connected branch which will be located on the minor 
diagonal itself. This area may not exist for small systems due to 
the very small number of such nodes. Note also that elements located 
to the left of the minor diagonal will appear high in the matrix 
structure in the case of small systems as opposed to the case of large 


systems. 


The shaded area is not full of off-diagonal elements. 
Rather, the concentration of these elements increases as one moves 
down along the main diagonal. This is because increasing row number 


almost corresponds to an increasing number of off-diagonal elements 
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An exception would be a case in which lightly connected nodes are 
connected together (e.g. chain circuit). In this case, heavily 
connected nodes will appear in the middle rows, resulting in a large 
number of fill in terms during the elimination process. The success 

of the method will, thus, depend on system structure, a characteristic 


of any practical ordering schemese 20, 


4.3 Storage Bookkeeping Technique 


Should system nodes be numbered as one unit, special 
programming techniques, apart from compact storage, would have not 
been necessary. System data, e.g. the bus admittance matrix, could 
be rearranged in the new order, and used directly. However, in the 
present method of solution, generator nodes and load nodes are 
handled separately, and, thus, should be so numbered. Voltage 
controlled busses have dual roles. They are control nodes since 
their voltage is a control variable. They are also dependent nodes, 
since their real powers are fixed, and the voltage angles are 
determined from the load-flow solution. Therefore, these nodes will 
participate in both node groups, and each of them will end up with 
two numbers to determine its order within each group. Rearranging 
the bus admittance network, say,to meet such a situation, would be 
impossible. 

To solve this problem, system data are left in the original 


arrangement. Correspondence arrays are formed. These arrays 


determine the original order of a given node, based on which group it 
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belongs to, and its order within that group. Once the original 
order of a node is determined, the Original data arrays can be 
readily used. Furthermore, the original node-branch connection 
array is split in two. The first, which is kept only for generator 
nodes*, gives, in the optimal order, the numbers of the generator 
nodes connected to each one. This array is used in computing and 
Storing the first and second order partial derivatives of the 
objective function (excluding the associated equality constraints) 


with respect to the control variables. 


The second consists of two separate arrays, one of which 
gives the list of the node-branch connections, in optimal order, 
within the group of load nodes, while the other gives the list of 
load nodes, in terms of their optimal order within their own group, 
connected to each generator node. The former array is used for 
computing and storing the Jacobian matrix, and the latter in 
computing and storing the first derivative of the cost function 
(excluding the associated equality constraints) with respect to the 
dependent variables. It is also used in computing the matrix [ag/ su] 
and its contribution in the Hessian matrix. 

To illustrate the just described technique, consider the 


IEEE 14-bus test system shown in Figure 4.2. Note that each of nodes 


ee ee ee ee ee ___ 


* It should be understood that reference to generator nodes, or 


load nodes, means the collection of such nodes plus voltage 


controlled nodes. 
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Figure 4.2 IEEE 14-Bus Test System Optimally Ordered 
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Figure 4.3 Storage Bookkeeping Arrangement 
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no. 4 and 5 (original order), which are voltage controlled nodes, 
has assumed two optimal orders. For instance, the node originally 
no. 5 is the second node in the group of generator nodes, while it 
is the third in the load nodes group. The correspondence arrays 
will, then, be given as in Figure 4.3 (a and b). Figures 4.3(a) 
gives the correspondence for generator nodes, while Figure 4.3(b) 
is for load nodes. In both arrays the position of an entry 
determines the node number in optimal order, while the entry is 


its number in the original order. 


The difference node-branch connection lists are shown in 
Figure 423 (c, di andie). Associated with each list 1s an index 
array which gives the starting position, in the respective connection 
list, of the list of nodes connected to each node of the node group 
concerned. For example, referring to Figure 4.3(c), which gives 
the node-branch connections within generator nodes, generator node 
no. 1 (optimal order; no. 3 in the original order) has no entry 
in the connection list in view of its isolation from other generator 
nodes. On the other hand, the list of generator nodes connected to 
generator node no. 5 (optimal order) starts at position 3 of the 
connection list. This means that node no. 5 is connected to nodes 


no. oP and: 4. 


4.4 Comparison of Ordering Schemes 


A comparison between the proposed method and schemes I 
and II of optimal ordering (see section 4.1) was carried out. The 


results of this comparison are presented in Tables 4.1 - 4.4. 
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Table 4.1 Ordering Times 


Required Storage 
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Scheme I 
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System Proposed 
Method 


30-bus 


57-bus 


118-bus 


Table 4.2 Storage of the Table of 


Factors 
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Table 4.1 gives the times required for ordering the 
three systems by the three methods. The times shown consist of the 
time required to order the two groups of nodes mentioned before, as 
well as all times required to load the program and read system data. 


These latter portions are the same for the three methods. 


It can be seen that scheme I is definitely the fastest 
of the three schemes tested. This is a well established resul thee 
While the proposed method and scheme II required about the same amount 
of time in the cases of the 30- and 118- bus systems, the proposed 
method was much faster in the 57- bus system case. This means 
that the speed of either the proposed scheme or scheme II jis not 
Only dependent on system size only, as in the case of scheme I, but 
also on system configuration. This should be expected since in 
both schemes ordering very much depends on how the system nodes are 


tied together. 


Table 4.2 gives the amount of storage required for the 
table of factors of the Jacobian matrix of the three systems using 
the three methods, as well as for the natural order. Storage required 
for the Hessian matrix is not shown since it represents a very smal] 


fraction of that required by the above mentioned table of factors*. 


a 


* In the case of natural order, which usually requires the largest 


amount of storage of all cases, the Hessian required only 26 


locations. 
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Storage Savings % 


Schemes | Proposed 
& Il Method 


30-bus 


57-bus 
118-bus 


Table 4.3 Percentage Savings in Storage 


Over that of Natural Order 


Scheme 
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Proposed 
Scheme 
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Table 4.4 Product of Ordering Time by 
Storage Required for Ordering 
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; og 
Storage requirement for the [ ae matrix is the same in all cases. 


Although the proposed method required between 14 and 27% 
more storage than that required by schemes I or II, the non-monotonic 
nature of such extra storage, backed up by the fact that scheme II 
was not better than scheme I in these cases, indicates that the 
required storage by an optimal scheme is dependent on system 


configuration, as mentioned before. 


Based on the amount of storage required by the natural 
order, the savings effected by the three schemes are shown in 
Table 4.3. These figures, coupled with the relative speed of each 
ordering scheme suggests that the proposed scheme iS,in general, 
comparable to scheme II, which is well accepted for load flow 


studies 20927 230534] 


The product of the time required for ordering by the 
amount of storage required is often used as a criterion to compare 
ordering penance ot This product is given in Table 4.4 for the 
three systems for scheme II and the proposed method. It backs up 
the above conclusion. 

A definite conclusion on the superiority of one scheme 
over the other should involve the consideration of a large number of 


power systems. This was not possible in the present Study. 


4.5 Minimum Cost Solution of the 30- Bus System 
The optimal solution of the minimum cost problem for the 


IEEE 30- bus test system is shown in Table 4.5. The system includes 
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3 generating nodes whose generation cost is the same as those of the 
14- bus system and are given in Table 3.1. Several starting 
conditions have been used for which the solution has been obtained in 
5 to 12 iterations. The 12 iterations figure refers to the flat 
start which, as mentioned before, is not a realistic operating point. 
This convergence rate confirms the earlier conclusion that it is 


independent of system size. 


Table 4.6 shows a comparison between the proposed solution 
method and the method of reference 22. Although the minimum 
generation cost as obtained by that latter method is lower, it will 
probably increase if better satisfaction of equality constraints 


is achieved. 
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Table 4.5 Minimum Cost Solution of the IEEE 30-Bus Test System 
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Cost of Generation = 1245.25 $/Hr. 


Table 4.6 Comparison Between the Proposed Method and the Method 


of Reference 22 for the 30-Bus System 


Minimum 
Cost 


5 1245.25 


Proposed <QeecxlUn 
Method 


3 1242.5 


Method of 1.5x10- 
Ref .22 


13. 15 


* Typical Mismatch: 0.3x10 ~ - 0.2x10— 
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CHAPTER V 


EFFECT OF PENALTY AND ACCELERATION FACTORS 


In this chapter, the effect of penalty and acceleration 
factors on the convergence of the proposed method is discussed. 
Since, practically, it is always possible, but not necessarily 
essential, to start the minimization process with a feasible operating 
point, only the effect of these two factors when a boundary of 


the feasible region is encountered from inside is considered. 


Furthermore, the concept of a fixed penalty factor 
developed with the method, is compared to the usual concept of a 


monotonically increasing sequence of penalty factors. 


Since the effects of penalty and acceleration factors 
complement one another, no attempt has been made to consider each 


separately. 


5.1 Effect of Penalty and Acceleration Factors on Convergence 


Referring to Figure 5.1, which represents schematically 


ax 


Figure 5.1 Effect of Penalty and Acceleration Factors 
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the minimization process, let curve A represent a boundary of the 
feasible region due some functional constraint. Let points a and a' 
represent, respectively, the unconstrained, and constrained optimal 
points. Also, let point b represent the Operating point at the 


Start of'a particular iteration. 


Assume that the unaccelerated move at that iteration 
moves the operating point b to a new point c outside the feasible 
region. If the penalty factor is small, point c is likely to have 
a lower cost (plus penalty) than point b, thus becoming the 
Starting point of the next iteration. Since the penalty factor is 
kept fixed, point c will continue to move toward point a, with very 
little attention to the satisfaction of functional constraints. 

No acceleration will likely be introduced, and eventually an un- 


feasible solution will be obtained. 


As long as the value of the penalty factor is not high 
enough to fully activate the acceleration mechanism, described in 
Chapter II, the minimization pattern will be similar to the one 
discussed above. Convergence, as measured by the number of iterations, 


will become worse as the penalty factor is increased. 


Suppose, now, that the value of the penalty factor is high 
enough, so that point c is of higher cost than point b. The acceleration 
mechanism will start functioning and the minimization move will be 
bd rather than bc. This pattern is continued, bd d'd"..., until the 
optimal point a' is located. In this case, constraint violation may 


be allowed, as long as the decrease in the objective function over- 
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comes such violation. Furthermore, a move from a point such as d" 
will try to produce a decrease in both the objective function and 
constraint violation. In such a way the Operating point will stay 


close to the boundary, resulting in a better convergence rate. 


If the penalty factor is higher than necessary, the 
minimization process may converge to a point such as e, which will, 
then, be a suboptimal solution. However, it will be characterized 
by a sudden or continued decrease in the value of the acceleration 
factor, indicating a premature convergence. The reason is that any 
Slight violation of a constraint will result in a large increase in 
the objective function, and the acceleration factor will continue to 
decrease to avoid such violation, and to keep the operating point 
within the feasible region. If the operating point at the start of 
the iteration is close to the boundary, it is likely that the value 
of the acceleration factor will drop to a very low level and could 
result in the satisfaction of the convergence criterion, without 
actually locating the optimal point. On the other hand, if constraint 
violation does occur in the early stages of the process, when the 
objective function changes rapidly, the emphasis will be on the 
satisfaction of the violated constraint, and very little attention is 
given to the minimization of the original objective function. The 
operating point will be pulled well inside the feasible region, far 


from the optimal solution. The result is a very slow convergence to 
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a suboptimal solution. This is shown by trajectory bby bo: in 


Bigure 5.1. 


This behaviour is demonstrated in Tables 5.1 and 5.2, which 
refer to two different cases of the 5- bus system. Similar behaviour 


is also experienced in the cases of the 14- bus system. 


Inelable 5.15 "1t-1s clear that, acceleration dic notetare 
full effect upto r=10", and thus, as expected convergence became 
worse as the penalty factor is increased. Constraint violation was 
smaller, however. In the case of r=10", acceleration took full 
action resulting in an excellent convergence behaviour in both rate 


and constraint violation. 


Similar behaviour is experienced in Table 5.2. However, 
in the case of rg? ns 10°, acceleration took effect during the iteration 
process rather than at its end. The fact that the process converged 
with unity acceleration factor and larger cost shows, as have also 
been discussed in Chapter III, that the cost contours are very flat 
and the satisfaction of a practical convergence criterion ,max su, <0", 
can occur in a region rather than a point. In this particular case a 


3 3 


range of values between 10” and 5xl0° for the penalty factor is 


satisfactory. 


The behaviour of the iteration process for rgclo" and 10° 


was typical of one with higher than necessary penalty factor. With 
the former value the operating point was repeatedly moved well inside 
the feasible region whenever constraint violation occur, whereas with 


the latter, it was not allowed to leave the feasible region in the 
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Table 5.1 Convergence Behaviour for Different Penalty Factors ; 


Min. Cost Problem, 5- Bus System, Free Gen. Voltages 


a NS ee eee ee 


Penalty No. of Aecr Violation of Cost at 

factor Iter. Factor Ue Constraint Solution 
= 2 7 -2 

ey aad 9 i 0 1.582x10 VS 

r=10° 10 1.0 1.351x10°¢ 757.22 

ry=10" Oscillatoryesoiution. Acc, Factoru=s055 

ry=10” 6 0.125 0.90x1074 757.78 


Table 5.2 Convergence Behaviour for Different Penalty Factors: 
Min. Cost Problem, Modified 5-Bus System, -0.4<Q,<0.4 


Fixed Gen. Voltages 


Penalty No. of Acc. Violation of COSteau 
Factor Iter. Factor Q3 Constraint Solution 
F700 7 1.0 7.30x107° 757.57 
rz 10° 7 1.0 1.03x10"2 757.63 
rgzl0” 9 1.0 Dasa” 75 7ics 
rg2x10™ 10 1.0 2.06x10"> 757H6T 
rg=5x10" 9 1.0 1.67x10"° 758.31 
rgrlo" 7 0.0625 1.72x107° 758.86 
rgn 10” 3 0.0625 : 759.42 
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first place. Both solutions are suboptimal. 


be FactorsAffecting the Choice of the Penalty Factor 


The main factor that determines whether or not acceleration 
will take place is the value of the penalty function due to violated 
constraints, and whether or not, when added to the objective function, 
it increases the value of such objective function over the previous 
value. Too large a penalty will certainly eat up any decrease that 
would be effected in the original objective function by the 
minimization process, forcing a premature convergence. On the other 
hand, too small a penalty will contribute a little, and a solution 


outside the feasible region would be obtained. 


Since the decrease in the objective function, particularly 
at the last stages of the process, is a very small fraction of that 
function, the penalty factor, which determines the value of the 
penalty function, should be chosen according to such a change. This 
means that it is dependent on the value of the objective function 
itself.. The larger this value is, the larger the penalty factor 


should be. 


Table 5.3 shows the convergence behaviour of two objective 
functions for varying penalty factor. Part (a) of the Table refers 
to the minimum cost problem of the 5- bus system when generator 
voltages are free. It is drawn from Table 5.1. Part (b) is obtained 
for the same case except the objective function was divided by a 
factor of 100. This factor should not affect the minimization process 


or the solution since it will cancel in the equation: 


. = 
a 
; “4 ai 
ce spin 1 fit ct 
i i I j ii 
_ ~4 - ’ 7 
i | ‘J if } ee wWiteiel G hi ’ i ) 3 4 aby 
i re wag + 
0 
a ee ee twtaat 20 
; ) r 
: } : ; aol ie 4 
, IPMS SVEN eG $e4 ay ot rs285190 
a Ye oe" 19a 
bd * 
. Or =é ely 7 ‘e yi ae Si i Tt a ; 
4 ¥ ‘ — gull ih J? ie nite scott a2 OT 
| + ey: 
rig vi ie i AT) ar a 
ib Gith SOs Bet cu rr a ( 
' q 5 ei 
: ; . Die 1 nal 
\ 
| ie +. ne fog Ay 
; i. nge f ig fi 0 ae 
i j ri’ > iit 7 ite. F r ; if Th 
J y | ve { ' * 
} | io i ’ J f 6 ; 
‘ ; ' } 2 ‘Tt ah a fi 7% ' 
i i 
ne 


: _ - 
YF ‘7 AMAaAr - reer i ‘, 7 7 4¢ 

H iy Vif WV aa?? ‘ Sas r Ay rr a F ai 4 ~ ir Ga af a 
ert: Nd Gps SA Te 4%) rT! ‘0g tong: e “ my “sh 
| Ps: en? 
7 7 aF S ‘ me pis veo Pi ut, 2 *\ ad mn yl he a 
pore vane fg “es : e a: or - 
eT esha zi a \ 11. RP wid a ma } he E ee x 
; » ee eee: 


~ 2 i. _ 
A aimee 
c cr iid eae 


107 


UOL}OUNY BALZIaLqO ayy 


UO 40ZLDe4 ALL PUddg 40 JDUSPUsdaq E°S a|qel, 


-01X6 a2 lo 9 
UOLINLOS AA0ZeLLLISO 
-01%6 GcL°0 9 2-OLX SEL Oak OL 
Wop NLS “Morel L9s0 S| Se0tXzeci Ou 6 
2-O LX LSE“ L Gel OL 
2-01X28S "| O°L 6 
UOLZELOLA 407Ie4 "492] UOLZELOLA 40YIC4 "4dt] 4o}oe4 
*xXeW "99 $O ‘ON *XOW Oy Ai [euad 


uOLZOUNY aALYDAaLaQ pajeds UOLLOUNY SALZOALAQ | eWwUON 


a 
an ena 


Yb 
fonu7 ev Fis JF bd bel: 
— I _ - 


> i 
Bi J aes 


108 


au = -H! vf (5.1) 


‘<7 tok 


which gives the minimization move. 


It can be seen that parts (a) and (b) would be identical 
if the factor of 100 is also taken into consideration in the value 
of the penalty factor. The same convergence behaviour is obtained 
for both objective functions when the penalty factor in both cases 


are related by this factor of 100. 


It should also be mentioned here that a penalty factor of 
200 was found acceptable for the minimum loss problem of the 5- bus 
system (Table 3.9 ) as compared to 1000 and 10° for the two cases 
presented here. The value of the objective function in the three 
cases is 1.644, 7.59 and 759.78 respectively. Correlation between 
the values of the peeerite function and the penalty factor is evident. 
The penalty factor is about 100 times larger than the objective function. 
This ratio was also satisfactory for all other cases for this type 


OF constraint. 


Another factor affecting the choice of the penalty factor 
is the type of constraint. It can be observed from Tables 5.1 and 
5.2 that while a value of 10° was acceptable for the voltage constraint 
at load node no. 5, a value of 10° tOmo” 10° was satisfactory for 
reactive generation constraint at voltage controlled node no. 3. 


This difference can be explained in the following way. 


One should first recognize that the former type of constraint 
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is a linear constraint of only one dependent variable, i.e. the 
voltage at a load node, where as the latter is a nonlinear constraint 
which is a function of both control and dependent variables. These 
are the two types of functional constraints involved in the optimal 


load flow problem. In view of the equations 


3(f+w heel 
ue ee d= 0 (5.2) 
aay i 
a( Ft 
vf set et ese (5.3) 


which represent the gradients of the penalized objective function 
with respect to dependent and control variables respectively*, it can 
be seen that the first type of constraint, which has only a first 
order derivative with respect to a particular dependent variable, 
will have only an indirect contribution, through the vector i, to 
the gradient vf of equation (5.3). On the other hand, the second 
type of pone cra ints will have both such indirect contribution and 

a direct one through the derivative ~. A lower penalty factor 

for that latter type should, then, is CREEL An acceptable value 
for the penalty factor in this case is about 2 to 10 times the value 
of the objective function. 

NII eee ee ee 


* These equations are identical to equations (2.18) and (2.19) 


except a penalty function w is added to the objective function 
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5.3 Comparison of Penalty Factor Schemes 


A proper mathematical way to solve a constrained 
minimization problem using penalty function approach is to start the 
minimization process with a small penalty factor and then increase 
its value after convergence is achieved. This process is continued 
until an acceptable solution is obtained. This process will, then, 
require a large number of iterations as it consists of several 
minimization subproblems. This is in contrast with only one such 
minimization problem, when the concept of a fixed penalty factor 


is used. 


If the objective function has a distinctive constrained 
optimal solution, it could be argued that the former method will be 
definitely the superior. However, this is not the case in minimization 
problems of power systems. The contours of the objective function 
are very flat, and any practical convergence criterion will locate a 
neighbourhood of the optimal solution rather than the solution point 
itself. Within such neighbourhood cost variation is minimal and, 
from a practical point of view, the extra effort required to locate 
the true optimum is not worthwhile. 

Consider Table 5.4 which shows the number of iterations 
required to obtain the solution of the minimum cost problem for the 
5- bus and the modified 5- bus systems using both penalty factor approaches. 
In the first case, generator voltages were free, while in the second 


case they were fixed and -0.4<0,<0.4. 
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Costeat 
Solution 


Fixed Penalty 


Did not converge in 35 iterations 


9 Deal = REGGE 
14 1.12x107°9 757.88 
119x100" 759.16 


Table 5.4 Comparison of Penalty Factor 


Fixed Penalty 


Factor 


Modified 5-Bus 


Approaches 


1 


7s 
ee 
J 


It is quite clear that.for all practical reasons, 
the concept of fixed penalty factor is more acceptable. In the 
first case the difference in optimal cost indeed does not justify 
the extra iterations. In the second case a varying penalty factor 
approach was not acceptable with regard to both convergence rate 


and the solution obtained. It does point out, however, the 
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sensitivity and cruciality of the growth rate of the penalty factor. 


aely nial 5 
| (0G “pen aa: 
Hohl . ony. Lem Ses A500 “vii : 83a 224 
| a ennas ot} al eee sed tt ‘i 
is ‘ye ai apo v. 


fiitty Porn ie ono’ >it 
ia 


CHAPTER VI 


CONCLUSIONS 


In this thesis the problem of optimal load flow of power 
systems is investigated. A brief review of the development of the 
problem and its solution methods is given. The Carpentier 


formulation of the exact economic dispatch problem, which is the 


basic starting point of all recent developments, is also presented. 


The solution method proposed in this thesis, which is 


also based on Carpentier formulation, is characterized by: 


1) All generators are treated as swing busses. 

2) Newton's method, using the Hessian matrix of the objective 
function, is used to compute the adjustments to the control 
variables. 


3) A simple acceleration scheme is used. 


4) A fixed penalty factor is used throughout the iterative process. 


By treating all generators as swing busses, several 


advantages have been achieved: 
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a) The sensitivity of some problems to the location of one particular 


swing bus is eliminated. 

b) Storage requirement is greatly reduced, as compared to other 
methods, due to the elimination of generator equations from the 
set of load flow equations. The Jacobian matrix of load flow 


equations is much smaller, resulting also in a large reduction 
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in computation efforts required to solve the load flow DORLION OT 
the problem. 

c) In addition, generator voltage magnitudes are kept free to adjust 
themselves within the prescribed range, to achieve true 
minimization. Although this is possible with other methods, 
the only results that have been reported are those of Billinton 


[20] 


and Sachdeva using a suboptimum technique of sequential 


real/reactive power dispatch. 


The use of Newton's method to adjust the set of control 
variables, as defined in this thesis, provided excellent convergence 
rate. The Hessian matrix needed for this purpose is of a low order 
and is extremely sparse, thus posing no problems in its handling or 


storage. 


By using the simple acceleration scheme described in this 
thesis, the complexities and efforts needed to evaluate an optimal 
acceleration factor is eliminated. Yet both may produce the same 


results. 


The fixed penalty factor approach, complimented by the 
acceleration scheme, proved to be superior to the usual way in which 
penalty factors are used. The process converges in fewer iterations 
with little or no sacrifice in the final cost. This latter point is 
practically of no concern because of the flatness of the cost contours 
and the impossibility of implementing, exactly, the true mathematical 
optimal solution. By relating the choice of this fixed penalty 


factor to the objective function and the type of constraint, a rough 
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guide to estimate its value was possible. 


The method is extensively tested using two standard test 
systems. Four basic problems are solved: the minimum cost 
problem, or the exact economic dispatch, the minimum loss problem, 
the minimum fuel problem, and the combined cost-fuel minimization 
problem. The first two problems are well established, whereas the 


last two are newly defined and formulated for power systems. 


An optimal ordering scheme, for use with large systems, 
is also developed in this thesis. Near optimal order is achieved 
by moving off diagonal elements to the right of the minor diagonal 
and as far below it as possible. The scheme was found to be 


comparable with two of the well accepted schemes. 


For future considerations, it is worthwhile investigating 
the feasibility of applying the method to the problem of hydro- 
thermal dispatch, in which optimization should be carried out over a 
time period, rather than at a particular instant. The method can 
also be investigated from the point of view of system planning and 
expansion. More work also need to be done on the fixed penalty 


factor approach to establish a firm guideline for the choice of the 


proper value. 
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Appendix A 


Given here are the elements of the different matrices and 


vectors involved in the solution method presented in Chapter II. 


A.1 The Jacobian Matrix 


Load flow equations at a load node are given by: 
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The elements of the Jacobian are the derivatives of these 
quantities with respect to the dependent variables x. Self 


derivatives are given by: 
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Similarly 
oq. 
] = = -s 
25, = C, V ee cos(6,.) (A.5) 
oq. 
<a 
Sie ey ee Valeep cuba (A.6) 


Mutual derivatives are given by: 


Op. 

‘if ; 
3, = eg ee aay (A.7) 
op. P 
wv, = Vita cos(6;-35-8;.) (A.8) 
oq. ; 
355 == ENE GSES OO (A.9) 
24, 
av, = Vita sin(s;-65-6;,) (A.10) - 


In equations (A.3) - (A.10) i and j stand for load busses. 
A.2 The Cost Function and its derivatives 


It is assumed that the cost function is a quadratic 


function in generated real powers as follows: 


GS ny ae as 8? tee de (A.11) 
feo ea eee 


id re Tee 
oy | eo iee autem | 
. 7 fa 


1) iat 


(sae be shit 


: ; 
‘ ce 
4 . 
} ruer 3 
. F ‘ | ftw - 
Pr! I + 
f ; 7 
a : ‘ ] 
. eer ys 
*E w eat ore ot | 
7) 
lt 
l ; 
t} TT ‘> e 
et ae 
ag 7 ; 
; One tae i 
od fon! wo) biwte |, hae Ot A) ~ fe een 
a 4 . 
: ; 
. mn 7 - _ 
Vt an yi} =u fs 90.6 a 


: aa arr we Tr as 4 
SVinti Lie ¢ a! 0 ‘tA. foe: nit Ynits fot we: a. 


“eae es dat i Pe 


where i is for a generator bus, a.» b. and c. 
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oe the real power generated at bus i, is given by: 
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A.2.1 Elements of the derivative vector os 
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A.2.2 Elements of the vector 2* 
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oP 
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t(D dae caeem) sin(s, - 6) - O1)] (A.22) 


In Equations (A.15) - (A.22), k and 2 are for generator 


busses, and the derivatives of generated powers a , with respect to 
j 
generator variables Sas Me % 8; and Ve are given by: 
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A.3 The Matrix [ag/au] 


Equations (A.7) - (A.10) give the elements of the matrix 
[ag/su] provided that subscript j is for a generator node rather than 
a load node. 
A.4 The Hessian Matrix 


The elements of the Hessian matrix are the derivatives of 


relations (2.19) with respect to the control variables. Bearing in 
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mind equation (A.7) - (A.10) and the provision of section AN 3), 


these elements reduce to: 
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where, in (A.27) - (A.32), j is for a load bus, and k and & are for, 
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A.5 Derivatives of the Penalty Function 


The penalty function (2.33) can be written, in more detail, 
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APPENDIX B 


DESCRIPTION OF TEST SYSTEMS 


Given below are the line diagrams and the bus admittance 
matrices, in polar form, of the systems used in this thesis. Also 
given are their loading conditions. In calculating the bus admittance 
matrices, off-nominal transformer tap ratios, and static capacitors 
have been taken into account. Regulated busses (generator and voltage 
controlled busses) in each system are indicated in the solutions given 


in Chapters 3 and 4. 
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Table B.2 


SYSTEM LOADS 
(IN P.U. ON A 100 MVA Base) 


REAL REACTIVE BUS REAL 
0.0 0.0 4 0.400 
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Figure B.1 Line Diagram of the 5-bus Test System 
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B.2 The IEEE 14-Bus Test System 


©) SYNCHRONOUS CAPACITORS 

© GENERATORS 

(a) BUS-CODE DIAGRAM 
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Figure B.2 Line Diagram of the IEEE 14- Bus Test 
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Table B.4 SYSTEMS SEORNDS 
(IN PLU. ON A 100 MVA BASE) 


BUS REAL REACTIVE BUS 
1 0.0 0.0 8 
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3 Oa uh2 OOS 10 
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Table B.5 DIAGONAL AND NON-ZERO UPPER DIAGONAL 
ELEMENTS OF THE BUS ADMITTANCE MATRIX 
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Table BiG SYSTEM LOADS 
(IN PRIUS SON TA 9100 NVA BASE) 


BUS REAL REACTIVE BUS REAL REACTIVE 
1 0.0 0.0 16 0.035 0.018 
2 One 7 Oeatze 17 0.090 0.058 
3 0.024 0.012 18 OS OS)” 0.009 
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B.4 The IEEE 57-Bus Test System 
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Figure B.4 Line Diagram of the IEEE 57-Bus Test System . 
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Table B.7 DIAGONAL AND NON-ZERO UPPER DIAGONAL 
ELEMENTS OF THE BUS ADMITTANCE MATRIX 
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Table= B28 GS YSTEM SAEOADS 
(IN P.U. ON A 100 MVA BASE) 


REAL REACTIVE BUS REAL REACTIVE 
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A line diagram of the IEEE 118-Bus test system can be 


found apenas ad, Only the bus admittance matrix and loads of the 


system are given here. 
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